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

Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development

This commit is contained in:
Martin Diehl 2016-07-04 21:43:33 +02:00
parent cd5eddb444 a317636e8b
commit 3ada327d3f
39 changed files with 749 additions and 12629 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
DAMASK_env.bat Normal file
View File

@ -0,0 +1,19 @@
:: sets up an environment for DAMASK on Windows
:: usage: call DAMASK_env.bat
@echo off
set LOCATION=%~dp0
set DAMASK_ROOT=%LOCATION%\DAMASK
set DAMASK_NUM_THREADS=2
chcp 1252
Title Düsseldorf Advanced Materials Simulation Kit - DAMASK, MPIE Düsseldorf
echo.
echo Düsseldorf Advanced Materials Simulation Kit - DAMASK
echo Max-Planck-Institut für Eisenforschung, Düsseldorf
echo http://damask.mpie.de
echo.
echo Preparing environment ...
echo DAMASK_ROOT=%DAMASK_ROOT%
echo DAMASK_NUM_THREADS=%DAMASK_NUM_THREADS%
set DAMASK_BIN=%DAMASK_ROOT%\bin
set PATH=%PATH%;%DAMASK_BIN%
set PYTHONPATH=%PYTHONPATH%;%DAMASK_ROOT%\lib

2
Makefile
View File

@ -21,11 +21,11 @@ marc:
processing:
@if hash cython 2>/dev/null; then \
cd ./lib/damask; \
ln -s orientation.py corientation.pyx; \
CC=gcc python setup_corientation.py build_ext --inplace; \
rm -rv build; \
rm *.c; \
fi
@./installation/compile_CoreModule.py ${MAKEFLAGS}
.PHONY: tidy
tidy:

2
VERSION
View File

@ -1 +1 @@
v2.0.0-297-ga27aba1
v2.0.0-341-gaf4307e

31
code/DAMASK_spectral.f90
View File

@ -617,12 +617,12 @@ program DAMASK_spectral
timeinc = timeinc/2.0_pReal
elseif (solres(1)%termIll) then ! material point model cannot find a solution, exit in any casy
call IO_warning(850_pInt)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written (e.g. for regridding)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
elseif (continueCalculation == 1_pInt) then
guess = .true. ! accept non converged BVP solution
else ! default behavior, exit if spectral solver does not converge
call IO_warning(850_pInt)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written (e.g. for regridding)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
endif
else
guess = .true. ! start guessing after first converged (sub)inc
@ -722,34 +722,29 @@ end program DAMASK_spectral
!> @brief quit subroutine to mimic behavior of FEM solvers
!> @details exits the Spectral solver and reports time and duration. Exit code 0 signals
!> everything went fine. Exit code 1 signals an error, message according to IO_error. Exit code
!> 2 signals request for regridding, increment of last saved restart information is written to
!> 2 signals no converged solution and increment of last saved restart information is written to
!> stderr. Exit code 3 signals no severe problems, but some increments did not converge
!--------------------------------------------------------------------------------------------------
subroutine quit(stop_id)
use prec, only: &
pInt
use numerics, only: &
worldrank
implicit none
integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
if (worldrank == 0_pInt) then
call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
endif
call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
if (stop_id == 0_pInt) stop 0 ! normal termination
if (stop_id < 0_pInt) then ! trigger regridding
if (worldrank == 0_pInt) &
write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)
if (stop_id < 0_pInt) then ! terminally ill, restart might help
write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)
stop 2
endif
if (stop_id == 3_pInt) stop 3 ! not all incs converged

10
code/FEsolving.f90
View File

@ -72,11 +72,9 @@ subroutine FE_init
integer(pInt), allocatable, dimension(:) :: chunkPos
#endif
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- FEsolving init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- FEsolving init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
modelName = getSolverJobName()
#ifdef Spectral
@ -153,10 +151,6 @@ subroutine FE_init
200 close(FILEUNIT)
endif
!--------------------------------------------------------------------------------------------------
! the following array are allocated by mesh.f90 and need to be deallocated in case of regridding
if (allocated(calcMode)) deallocate(calcMode)
if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
#endif
if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0_pInt) then
write(6,'(a21,l1)') ' restart writing: ', restartWrite

21
code/IO.f90
View File

@ -80,25 +80,10 @@ subroutine IO_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
implicit none
integer(pInt) :: worldrank = 0_pInt
#ifdef PETSc
#include <petsc/finclude/petscsys.h>
PetscErrorCode :: ierr
#endif
external :: &
MPI_Comm_rank, &
MPI_Abort
#ifdef PETSc
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
#endif
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- IO init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- IO init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
end subroutine IO_init
@ -1640,8 +1625,6 @@ subroutine IO_error(error_ID,el,ip,g,ext_msg)
msg = 'update of gamma operator not possible when pre-calculated'
case (880_pInt)
msg = 'mismatch of microstructure count and a*b*c in geom file'
case (890_pInt)
msg = 'invalid input for regridding'
case (891_pInt)
msg = 'unknown solver type selected'
case (892_pInt)

7
code/Makefile
View File

@ -351,7 +351,7 @@ DAMASK_spectral.o: INTERFACENAME := spectral_interface.f90
SPECTRAL_SOLVER_FILES = spectral_mech_AL.o spectral_mech_Basic.o spectral_mech_Polarisation.o \
spectral_thermal.o spectral_damage.o
SPECTRAL_FILES = C_routines.o system_routines.o prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o \
SPECTRAL_FILES = C_routines.o system_routines.o prec.o DAMASK_interface.o IO.o numerics.o debug.o math.o \
FEsolving.o mesh.o material.o lattice.o \
$(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \
crystallite.o \
@ -401,7 +401,7 @@ DAMASK_FEM.exe: INCLUDE_DIRS += -I./
FEM_SOLVER_FILES = FEM_mech.o FEM_thermal.o FEM_damage.o FEM_vacancyflux.o FEM_porosity.o FEM_hydrogenflux.o
FEM_FILES = prec.o DAMASK_interface.o FEZoo.o IO.o libs.o numerics.o debug.o math.o \
FEM_FILES = prec.o DAMASK_interface.o FEZoo.o IO.o numerics.o debug.o math.o \
FEsolving.o mesh.o material.o lattice.o \
$(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \
crystallite.o \
@ -612,9 +612,6 @@ debug.o: debug.f90 \
numerics.o
numerics.o: numerics.f90 \
libs.o
libs.o: libs.f90 \
IO.o
IO.o: IO.f90 \

1
code/commercialFEM_fileList.f90
View File

@ -4,7 +4,6 @@
!> @details List of files needed by MSC.Marc, Abaqus/Explicit, and Abaqus/Standard
!--------------------------------------------------------------------------------------------------
#include "IO.f90"
#include "libs.f90"
#include "numerics.f90"
#include "debug.f90"
#include "math.f90"

12
code/core_quit.f90
View File

@ -1,12 +0,0 @@
!********************************************************************
! quit subroutine to satisfy IO_error for core module
!
!********************************************************************
subroutine quit(stop_id)
use prec, only: &
pInt
implicit none
integer(pInt), intent(in) :: stop_id
end subroutine

12
code/libs.f90
View File

@ -1,12 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief dummy source for inclusion of Library files
!--------------------------------------------------------------------------------------------------
module libs
!nothing in here
end module libs
#include "../lib/IR_Precision.f90"
#include "../lib/Lib_Base64.f90"
#include "../lib/Lib_VTK_IO.f90"

19
code/math.f90
View File

@ -72,10 +72,6 @@ module math
3_pInt,3_pInt &
],[2,9]) !< arrangement in Plain notation
#ifdef Spectral
include 'fftw3.f03'
#endif
public :: &
math_init, &
math_qsort, &
@ -163,21 +159,6 @@ module math
math_rotate_forward33, &
math_rotate_backward33, &
math_rotate_forward3333
#ifdef Spectral
public :: &
fftw_set_timelimit, &
fftw_plan_dft_3d, &
fftw_plan_many_dft_r2c, &
fftw_plan_many_dft_c2r, &
fftw_plan_with_nthreads, &
fftw_init_threads, &
fftw_alloc_complex, &
fftw_execute_dft, &
fftw_execute_dft_r2c, &
fftw_execute_dft_c2r, &
fftw_destroy_plan, &
math_tensorAvg
#endif
private :: &
math_partition, &
halton, &

597
code/mesh.f90
View File

@ -116,12 +116,8 @@ module mesh
#endif
#ifdef Spectral
#ifdef PETSc
#include <petsc/finclude/petscsys.h>
include 'fftw3-mpi.f03'
#else
include 'fftw3.f03'
#endif
#endif
! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS)
@ -413,18 +409,13 @@ module mesh
mesh_build_ipVolumes, &
mesh_build_ipCoordinates, &
mesh_cellCenterCoordinates, &
mesh_init_postprocessing, &
mesh_get_Ncellnodes, &
mesh_get_unitlength, &
mesh_get_nodeAtIP
#ifdef Spectral
public :: &
mesh_spectral_getGrid, &
mesh_spectral_getSize, &
mesh_nodesAroundCentres, &
mesh_deformedCoordsFFT, &
mesh_volumeMismatch, &
mesh_shapeMismatch
mesh_spectral_getSize
#endif
private :: &
@ -436,8 +427,7 @@ module mesh
mesh_spectral_build_nodes, &
mesh_spectral_build_elements, &
mesh_spectral_build_ipNeighborhood, &
#endif
#ifdef Marc4DAMASK
#elif defined Marc4DAMASK
mesh_marc_get_tableStyles, &
mesh_marc_count_nodesAndElements, &
mesh_marc_count_elementSets, &
@ -448,8 +438,7 @@ module mesh
mesh_marc_build_nodes, &
mesh_marc_count_cpSizes, &
mesh_marc_build_elements, &
#endif
#ifdef Abaqus
#elif defined Abaqus
mesh_abaqus_count_nodesAndElements, &
mesh_abaqus_count_elementSets, &
mesh_abaqus_count_materials, &
@ -473,11 +462,7 @@ module mesh
mesh_tell_statistics, &
FE_mapElemtype, &
mesh_faceMatch, &
mesh_build_FEdata, &
mesh_write_cellGeom, &
mesh_write_elemGeom, &
mesh_write_meshfile, &
mesh_read_meshfile
mesh_build_FEdata
contains
@ -487,9 +472,7 @@ contains
!! Order and routines strongly depend on type of solver
!--------------------------------------------------------------------------------------------------
subroutine mesh_init(ip,el)
#ifdef Spectral
use, intrinsic :: iso_c_binding
#endif
use DAMASK_interface
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use IO, only: &
@ -531,11 +514,9 @@ subroutine mesh_init(ip,el)
integer(pInt) :: j
logical :: myDebug
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem)
if (allocated(mesh_mapFEtoCPnode)) deallocate(mesh_mapFEtoCPnode)
@ -562,25 +543,18 @@ subroutine mesh_init(ip,el)
myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt)
#ifdef Spectral
#ifdef PETSc
call fftw_mpi_init()
#endif
call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file...
if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6)
grid = mesh_spectral_getGrid(fileUnit)
geomSize = mesh_spectral_getSize(fileUnit)
#ifdef PETSc
gridMPI = int(grid,C_INTPTR_T)
alloc_local = fftw_mpi_local_size_3d(gridMPI(3), gridMPI(2), gridMPI(1)/2 +1, &
MPI_COMM_WORLD, local_K, local_K_offset)
grid3 = int(local_K,pInt)
grid3Offset = int(local_K_offset,pInt)
size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal)
size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal)
#endif
if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6)
call mesh_spectral_count()
if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6)
@ -592,8 +566,7 @@ subroutine mesh_init(ip,el)
if (myDebug) write(6,'(a)') ' Built nodes'; flush(6)
call mesh_spectral_build_elements(FILEUNIT)
if (myDebug) write(6,'(a)') ' Built elements'; flush(6)
#endif
#ifdef Marc4DAMASK
#elif defined Marc4DAMASK
call IO_open_inputFile(FILEUNIT,modelName) ! parse info from input file...
if (myDebug) write(6,'(a)') ' Opened input file'; flush(6)
call mesh_marc_get_tableStyles(FILEUNIT)
@ -616,8 +589,7 @@ subroutine mesh_init(ip,el)
if (myDebug) write(6,'(a)') ' Counted CP sizes'; flush(6)
call mesh_marc_build_elements(FILEUNIT)
if (myDebug) write(6,'(a)') ' Built elements'; flush(6)
#endif
#ifdef Abaqus
#elif defined Abaqus
call IO_open_inputFile(FILEUNIT,modelName) ! parse info from input file...
if (myDebug) write(6,'(a)') ' Opened input file'; flush(6)
noPart = IO_abaqus_hasNoPart(FILEUNIT)
@ -666,15 +638,12 @@ subroutine mesh_init(ip,el)
if (myDebug) write(6,'(a)') ' Built shared elements'; flush(6)
call mesh_build_ipNeighborhood
#else
call mesh_spectral_build_ipNeighborhood(FILEUNIT)
call mesh_spectral_build_ipNeighborhood
#endif
if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6)
if (worldrank == 0_pInt) then
call mesh_tell_statistics
call mesh_write_meshfile
call mesh_write_cellGeom
call mesh_write_elemGeom
endif
if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) &
@ -963,7 +932,7 @@ subroutine mesh_build_ipCoordinates
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
myCoords = myCoords + mesh_cellnode(1:3,mesh_cell(n,i,e))
enddo
mesh_ipCoordinates(1:3,i,e) = myCoords / real(FE_NcellnodesPerCell(c),pReal)
mesh_ipCoordinates(1:3,i,e) = myCoords / FE_NcellnodesPerCell(c)
enddo
enddo
!$OMP END PARALLEL DO
@ -990,7 +959,7 @@ pure function mesh_cellCenterCoordinates(ip,el)
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el))
enddo
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / real(FE_NcellnodesPerCell(c),pReal)
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / FE_NcellnodesPerCell(c)
end function mesh_cellCenterCoordinates
@ -1417,11 +1386,9 @@ end subroutine mesh_spectral_build_elements
!> @brief build neighborhood relations for spectral
!> @details assign globals: mesh_ipNeighborhood
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_build_ipNeighborhood(fileUnit)
subroutine mesh_spectral_build_ipNeighborhood
implicit none
integer(pInt), intent(in) :: &
fileUnit
integer(pInt) :: &
x,y,z, &
e
@ -1558,332 +1525,8 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
nodes = nodes/8.0_pReal
end function mesh_nodesAroundCentres
!--------------------------------------------------------------------------------------------------
!> @brief calculate coordinates in current configuration for given defgrad
! using integration in Fourier space
!--------------------------------------------------------------------------------------------------
function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords)
use IO, only: &
IO_error
use numerics, only: &
fftw_timelimit, &
fftw_planner_flag
use debug, only: &
debug_mesh, &
debug_level, &
debug_levelBasic
use math, only: &
PI, &
math_mul33x3
implicit none
real(pReal), intent(in), dimension(:,:,:,:,:) :: F
real(pReal), dimension(3,size(F,3),size(F,4),size(F,5)) :: coords
real(pReal), intent(in), dimension(3) :: gDim
real(pReal), intent(in), dimension(3,3), optional :: FavgIn
real(pReal), intent(in), dimension(3), optional :: scalingIn
! allocatable arrays for fftw c routines
type(C_PTR) :: planForth, planBack
type(C_PTR) :: coords_fftw, defgrad_fftw
real(pReal), dimension(:,:,:,:,:), pointer :: F_real
complex(pReal), dimension(:,:,:,:,:), pointer :: F_fourier
real(pReal), dimension(:,:,:,:), pointer :: coords_real
complex(pReal), dimension(:,:,:,:), pointer :: coords_fourier
! other variables
integer(pInt) :: i, j, k, m, res1Red
integer(pInt), dimension(3) :: k_s, iRes
real(pReal), dimension(3) :: scaling, step, offset_coords, integrator
real(pReal), dimension(3,3) :: Favg
integer(pInt), dimension(2:3,2) :: Nyquist ! highest frequencies to be removed (1 if even, 2 if odd)
if (present(scalingIn)) then
where (scalingIn < 0.0_pReal) ! invalid values. in case of f2py -1 if not present
scaling = [1.0_pReal,1.0_pReal,1.0_pReal]
elsewhere
scaling = scalingIn
end where
else
scaling = 1.0_pReal
endif
iRes = [size(F,3),size(F,4),size(F,5)]
integrator = gDim / 2.0_pReal / PI ! see notes where it is used
res1Red = iRes(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
step = gDim/real(iRes, pReal)
Nyquist(2,1:2) = [iRes(2)/2_pInt + 1_pInt, iRes(2)/2_pInt + 1_pInt + mod(iRes(2),2_pInt)]
Nyquist(3,1:2) = [iRes(3)/2_pInt + 1_pInt, iRes(3)/2_pInt + 1_pInt + mod(iRes(3),2_pInt)]
!--------------------------------------------------------------------------------------------------
! report
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Restore geometry using FFT-based integration'
write(6,'(a,3(i12 ))') ' grid a b c: ', iRes
write(6,'(a,3(es12.5))') ' size x y z: ', gDim
endif
!--------------------------------------------------------------------------------------------------
! sanity check
if (pReal /= C_DOUBLE .or. pInt /= C_INT) &
call IO_error(0_pInt,ext_msg='Fortran to C in mesh_deformedCoordsFFT')
!--------------------------------------------------------------------------------------------------
! allocation and FFTW initialization
defgrad_fftw = fftw_alloc_complex(int(res1Red *iRes(2)*iRes(3)*9_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8)
coords_fftw = fftw_alloc_complex(int(res1Red *iRes(2)*iRes(3)*3_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8)
call c_f_pointer(defgrad_fftw, F_real, &
[iRes(1)+2_pInt-mod(iRes(1),2_pInt),iRes(2),iRes(3),3_pInt,3_pInt])
call c_f_pointer(defgrad_fftw, F_fourier, &
[res1Red, iRes(2),iRes(3),3_pInt,3_pInt])
call c_f_pointer(coords_fftw, coords_real, &
[iRes(1)+2_pInt-mod(iRes(1),2_pInt),iRes(2),iRes(3),3_pInt])
call c_f_pointer(coords_fftw, coords_fourier, &
[res1Red, iRes(2),iRes(3),3_pInt])
call fftw_set_timelimit(fftw_timelimit)
planForth = fftw_plan_many_dft_r2c(3_pInt,[iRes(3),iRes(2) ,iRes(1)],9_pInt,& ! dimensions , length in each dimension in reversed order
F_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt-mod(iRes(1),2_pInt)],& ! input data , physical length in each dimension in reversed order
1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt-mod(iRes(1),2_pInt)),& ! striding , product of physical lenght in the 3 dimensions
F_fourier,[iRes(3),iRes(2) ,res1Red],&
1_pInt, iRes(3)*iRes(2)* res1Red,fftw_planner_flag)
planBack = fftw_plan_many_dft_c2r(3_pInt,[iRes(3),iRes(2) ,iRes(1)],3_pInt,&
coords_fourier,[iRes(3),iRes(2) ,res1Red],&
1_pInt, iRes(3)*iRes(2)* res1Red,&
coords_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt-mod(iRes(1),2_pInt)],&
1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt-mod(iRes(1),2_pInt)),&
fftw_planner_flag)
F_real(1:iRes(1),1:iRes(2),1:iRes(3),1:3,1:3) = &
reshape(F,[iRes(1),iRes(2),iRes(3),3,3], order = [4,5,1,2,3]) ! F_real is overwritten during plan creatio, is larger (padding) and has different order
!--------------------------------------------------------------------------------------------------
! FFT
call fftw_execute_dft_r2c(planForth, F_real, F_fourier)
!--------------------------------------------------------------------------------------------------
! if no average F is given, compute it in Fourier space
if (present(FavgIn)) then
if (all(FavgIn < 0.0_pReal)) then
Favg = real(F_fourier(1,1,1,1:3,1:3),pReal)/real(product(iRes),pReal) !the f2py way to tell it is not present
else
Favg = FavgIn
endif
else
Favg = real(F_fourier(1,1,1,1:3,1:3),pReal)/real(product(iRes),pReal)
endif
!--------------------------------------------------------------------------------------------------
! remove highest frequency in each direction, in third direction only if not 2D
if(iRes(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
F_fourier (res1Red, 1:iRes(2), 1:iRes(3), 1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal)
if(iRes(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
F_fourier (1:res1Red,Nyquist(2,1):Nyquist(2,2),1:iRes(3), 1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal)
if(iRes(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
F_fourier (1:res1Red,1:iRes(2), Nyquist(3,1):Nyquist(3,2),1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal)
!--------------------------------------------------------------------------------------------------
! integration in Fourier space
coords_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
do k = 1_pInt, iRes(3)
k_s(3) = k-1_pInt
if(k > iRes(3)/2_pInt+1_pInt) k_s(3) = k_s(3)-iRes(3)
do j = 1_pInt, iRes(2)
k_s(2) = j-1_pInt
if(j > iRes(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-iRes(2)
do i = 1_pInt, res1Red
k_s(1) = i-1_pInt
do m = 1_pInt,3_pInt
coords_fourier(i,j,k,m) = sum(F_fourier(i,j,k,m,1:3)*&
cmplx(0.0_pReal,real(k_s,pReal)*integrator,pReal))
enddo
if (any(k_s /= 0_pInt)) coords_fourier(i,j,k,1:3) = &
coords_fourier(i,j,k,1:3) / cmplx(-sum(k_s*k_s),0.0_pReal,pReal)
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! iFFT and freeing memory
call fftw_execute_dft_c2r(planBack,coords_fourier,coords_real)
coords = reshape(coords_real(1:iRes(1),1:iRes(2),1:iRes(3),1:3), [3,iRes(1),iRes(2),iRes(3)], &
order = [2,3,4,1])/real(product(iRes),pReal) ! weight and change order
call fftw_destroy_plan(planForth)
call fftw_destroy_plan(planBack)
call fftw_free(defgrad_fftw)
call fftw_free(coords_fftw)
!--------------------------------------------------------------------------------------------------
! add average to scaled fluctuation and put (0,0,0) on (0,0,0)
offset_coords = math_mul33x3(F(1:3,1:3,1,1,1),step/2.0_pReal) - scaling*coords(1:3,1,1,1)
forall(k = 1_pInt:iRes(3), j = 1_pInt:iRes(2), i = 1_pInt:iRes(1)) &
coords(1:3,i,j,k) = scaling(1:3)*coords(1:3,i,j,k) &
+ offset_coords &
+ math_mul33x3(Favg,step*real([i,j,k]-1_pInt,pReal))
end function mesh_deformedCoordsFFT
!--------------------------------------------------------------------------------------------------
!> @brief calculates the mismatch between volume of reconstructed (compatible) cube and
! determinant of defgrad at the FP
!--------------------------------------------------------------------------------------------------
function mesh_volumeMismatch(gDim,F,nodes) result(vMismatch)
use IO, only: &
IO_error
use debug, only: &
debug_mesh, &
debug_level, &
debug_levelBasic
use math, only: &
math_det33, &
math_volTetrahedron
implicit none
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
F
real(pReal), dimension(size(F,3),size(F,4),size(F,5)) :: &
vMismatch
real(pReal), intent(in), dimension(:,:,:,:) :: &
nodes
real(pReal), dimension(3) :: &
gDim
integer(pInt), dimension(3) :: &
iRes
real(pReal), dimension(3,8) :: coords
integer(pInt) :: i,j,k
real(pReal) :: volInitial
iRes = [size(F,3),size(F,4),size(F,5)]
volInitial = product(gDim)/real(product(iRes), pReal)
!--------------------------------------------------------------------------------------------------
! report and check
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Calculating volume mismatch'
write(6,'(a,3(i12 ))') ' grid a b c: ', iRes
write(6,'(a,3(es12.5))') ' size x y z: ', gDim
endif
if (any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]))&
call IO_error(0_pInt,ext_msg='Arrays F and nodes in mesh_volumeMismatch')
!--------------------------------------------------------------------------------------------------
! calculate actual volume and volume resulting from deformation gradient
do k = 1_pInt,iRes(3)
do j = 1_pInt,iRes(2)
do i = 1_pInt,iRes(1)
coords(1:3,1) = nodes(1:3,i, j, k )
coords(1:3,2) = nodes(1:3,i+1_pInt,j, k )
coords(1:3,3) = nodes(1:3,i+1_pInt,j+1_pInt,k )
coords(1:3,4) = nodes(1:3,i, j+1_pInt,k )
coords(1:3,5) = nodes(1:3,i, j, k+1_pInt)
coords(1:3,6) = nodes(1:3,i+1_pInt,j, k+1_pInt)
coords(1:3,7) = nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt)
coords(1:3,8) = nodes(1:3,i, j+1_pInt,k+1_pInt)
vMismatch(i,j,k) = &
abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,8),coords(1:3,4))) &
+ abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,8),coords(1:3,5))) &
+ abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,3),coords(1:3,4))) &
+ abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,3),coords(1:3,2))) &
+ abs(math_volTetrahedron(coords(1:3,7),coords(1:3,5),coords(1:3,2),coords(1:3,6))) &
+ abs(math_volTetrahedron(coords(1:3,7),coords(1:3,5),coords(1:3,2),coords(1:3,1)))
vMismatch(i,j,k) = vMismatch(i,j,k)/math_det33(F(1:3,1:3,i,j,k))
enddo; enddo; enddo
vMismatch = vMismatch/volInitial
end function mesh_volumeMismatch
!--------------------------------------------------------------------------------------------------
!> @brief Routine to calculate the mismatch between the vectors from the central point to
! the corners of reconstructed (combatible) volume element and the vectors calculated by deforming
! the initial volume element with the current deformation gradient
!--------------------------------------------------------------------------------------------------
function mesh_shapeMismatch(gDim,F,nodes,centres) result(sMismatch)
use IO, only: &
IO_error
use debug, only: &
debug_mesh, &
debug_level, &
debug_levelBasic
use math, only: &
math_mul33x3
implicit none
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
F
real(pReal), dimension(size(F,3),size(F,4),size(F,5)) :: &
sMismatch
real(pReal), intent(in), dimension(:,:,:,:) :: &
nodes, &
centres
real(pReal), dimension(3) :: &
gDim, &
fRes
integer(pInt), dimension(3) :: &
iRes
real(pReal), dimension(3,8) :: coordsInitial
integer(pInt) i,j,k
iRes = [size(F,3),size(F,4),size(F,5)]
fRes = real(iRes,pReal)
!--------------------------------------------------------------------------------------------------
! report and check
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Calculating shape mismatch'
write(6,'(a,3(i12 ))') ' grid a b c: ', iRes
write(6,'(a,3(es12.5))') ' size x y z: ', gDim
endif
if(any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]) .or.&
any([iRes/=size(centres,2), iRes/=size(centres,3), iRes/=size(centres,4)]))&
call IO_error(0_pInt,ext_msg='Arrays F and nodes/centres in mesh_shapeMismatch')
!--------------------------------------------------------------------------------------------------
! initial positions
coordsInitial(1:3,1) = [-gDim(1)/fRes(1),-gDim(2)/fRes(2),-gDim(3)/fRes(3)]
coordsInitial(1:3,2) = [+gDim(1)/fRes(1),-gDim(2)/fRes(2),-gDim(3)/fRes(3)]
coordsInitial(1:3,3) = [+gDim(1)/fRes(1),+gDim(2)/fRes(2),-gDim(3)/fRes(3)]
coordsInitial(1:3,4) = [-gDim(1)/fRes(1),+gDim(2)/fRes(2),-gDim(3)/fRes(3)]
coordsInitial(1:3,5) = [-gDim(1)/fRes(1),-gDim(2)/fRes(2),+gDim(3)/fRes(3)]
coordsInitial(1:3,6) = [+gDim(1)/fRes(1),-gDim(2)/fRes(2),+gDim(3)/fRes(3)]
coordsInitial(1:3,7) = [+gDim(1)/fRes(1),+gDim(2)/fRes(2),+gDim(3)/fRes(3)]
coordsInitial(1:3,8) = [-gDim(1)/fRes(1),+gDim(2)/fRes(2),+gDim(3)/fRes(3)]
coordsInitial = coordsInitial/2.0_pReal
!--------------------------------------------------------------------------------------------------
! compare deformed original and deformed positions to actual positions
do k = 1_pInt,iRes(3)
do j = 1_pInt,iRes(2)
do i = 1_pInt,iRes(1)
sMismatch(i,j,k) = &
sqrt(sum((nodes(1:3,i, j, k ) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,1)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i+1_pInt,j, k ) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,2)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i+1_pInt,j+1_pInt,k ) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,3)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i, j+1_pInt,k ) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,4)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i, j, k+1_pInt) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,5)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i+1_pInt,j, k+1_pInt) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,6)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,7)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i, j+1_pInt,k+1_pInt) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,8)))**2.0_pReal))
enddo; enddo; enddo
end function mesh_shapeMismatch
#endif
#ifdef Marc4DAMASK
!--------------------------------------------------------------------------------------------------
!> @brief Figures out table styles (Marc only) and stores to 'initialcondTableStyle' and
@ -3070,6 +2713,7 @@ use IO, only: &
implicit none
integer(pInt), intent(in) :: fileUnit
#ifndef Spectral
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) chunk, Nchunks
@ -3081,9 +2725,10 @@ use IO, only: &
mesh_periodicSurface = .true.
#else
mesh_periodicSurface = .false.
#if defined(Marc4DAMASK)
#ifdef Marc4DAMASK
keyword = '$damask'
#elif defined(Abaqus)
#endif
#ifdef Abaqus
keyword = '**damask'
#endif
@ -3691,7 +3336,6 @@ integer(pInt) function FE_mapElemtype(what)
'c3d20t')
FE_mapElemtype = 13_pInt ! Three-dimensional Arbitrarily Distorted quadratic hexahedral
case default
FE_mapElemtype = -1_pInt ! error return
call IO_error(error_ID=190_pInt,ext_msg=IO_lc(what))
end select
@ -3700,7 +3344,6 @@ end function FE_mapElemtype
!--------------------------------------------------------------------------------------------------
!> @brief find face-matching element of same type
!> @details currently not used, check if needed for HDF5 output, otherwise delete
!--------------------------------------------------------------------------------------------------
subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace)
@ -4511,212 +4154,6 @@ subroutine mesh_build_FEdata
end subroutine mesh_build_FEdata
!--------------------------------------------------------------------------------------------------
!> @brief writes out initial cell geometry
!--------------------------------------------------------------------------------------------------
subroutine mesh_write_cellGeom
use DAMASK_interface, only: &
getSolverJobName, &
getSolverWorkingDirectoryName
use IR_Precision, only: &
I4P
use Lib_VTK_IO, only: &
VTK_ini, &
VTK_geo, &
VTK_con, &
VTK_end
implicit none
integer(I4P), dimension(1:mesh_Ncells) :: celltype
integer(I4P), dimension(mesh_Ncells*(1_pInt+FE_maxNcellnodesPerCell)) :: cellconnection
integer(I4P):: error
integer(I4P):: g, c, e, CellID, i, j
cellID = 0_pInt
j = 0_pInt
do e = 1_pInt, mesh_NcpElems ! loop over cpElems
g = FE_geomtype(mesh_element(2_pInt,e)) ! get geometry type
c = FE_celltype(g) ! get cell type
do i = 1_pInt,FE_Nips(g) ! loop over ips=cells in this element
cellID = cellID + 1_pInt
celltype(cellID) = MESH_VTKCELLTYPE(c)
cellconnection(j+1_pInt:j+FE_NcellnodesPerCell(c)+1_pInt) &
= [FE_NcellnodesPerCell(c),mesh_cell(1:FE_NcellnodesPerCell(c),i,e)-1_pInt] ! number of cellnodes per cell & list of global cellnode IDs belnging to this cell (cellnode counting starts at 0)
j = j + FE_NcellnodesPerCell(c) + 1_pInt
enddo
enddo
error=VTK_ini(output_format = 'ASCII', &
title=trim(getSolverJobName())//' cell mesh', &
filename = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'_ipbased.vtk', &
mesh_topology = 'UNSTRUCTURED_GRID')
!ToDo: check error here
error=VTK_geo(NN = int(mesh_Ncellnodes,I4P), &
X = mesh_cellnode(1,1:mesh_Ncellnodes), &
Y = mesh_cellnode(2,1:mesh_Ncellnodes), &
Z = mesh_cellnode(3,1:mesh_Ncellnodes))
!ToDo: check error here
error=VTK_con(NC = int(mesh_Ncells,I4P), &
connect = cellconnection(1:j), &
!ToDo: check error here
cell_type = celltype)
error=VTK_end()
!ToDo: check error here
end subroutine mesh_write_cellGeom
!--------------------------------------------------------------------------------------------------
!> @brief writes out initial element geometry
!--------------------------------------------------------------------------------------------------
subroutine mesh_write_elemGeom
use DAMASK_interface, only: &
getSolverJobName, &
getSolverWorkingDirectoryName
use IR_Precision, only: &
I4P
use Lib_VTK_IO, only: &
VTK_ini, &
VTK_geo, &
VTK_con, &
VTK_end
implicit none
integer(I4P), dimension(1:mesh_NcpElems) :: elemtype
integer(I4P), dimension(mesh_NcpElems*(1_pInt+FE_maxNnodes)) :: elementconnection
integer(I4P):: error
integer(pInt):: e, t, n, i
i = 0_pInt
do e = 1_pInt, mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get element type
elemtype(e) = MESH_VTKELEMTYPE(t)
elementconnection(i+1_pInt) = FE_Nnodes(t) ! number of nodes per element
do n = 1_pInt,FE_Nnodes(t)
elementconnection(i+1_pInt+n) = mesh_element(4_pInt+n,e) - 1_pInt ! global node ID of node that belongs to this element (node counting starts at 0)
enddo
i = i + 1_pInt + FE_Nnodes(t)
enddo
error=VTK_ini(output_format = 'ASCII', &
title=trim(getSolverJobName())//' element mesh', &
filename = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'_nodebased.vtk', &
mesh_topology = 'UNSTRUCTURED_GRID')
!ToDo: check error here
error=VTK_geo(NN = int(mesh_Nnodes,I4P), &
X = mesh_node0(1,1:mesh_Nnodes), &
Y = mesh_node0(2,1:mesh_Nnodes), &
Z = mesh_node0(3,1:mesh_Nnodes))
!ToDo: check error here
error=VTK_con(NC = int(mesh_Nelems,I4P), &
connect = elementconnection(1:i), &
cell_type = elemtype)
!ToDo: check error here
error =VTK_end()
!ToDo: check error here
end subroutine mesh_write_elemGeom
!--------------------------------------------------------------------------------------------------
!> @brief writes description file for mesh
!--------------------------------------------------------------------------------------------------
subroutine mesh_write_meshfile
use IO, only: &
IO_write_jobFile
implicit none
integer(pInt), parameter :: fileUnit = 223_pInt
integer(pInt) :: e,i,t,g,c,n
call IO_write_jobFile(fileUnit,'mesh')
write(fileUnit,'(A16,E10.3)') 'unitlength', mesh_unitlength
write(fileUnit,'(A16,I10)') 'maxNcellnodes', mesh_maxNcellnodes
write(fileUnit,'(A16,I10)') 'maxNips', mesh_maxNips
write(fileUnit,'(A16,I10)') 'maxNnodes', mesh_maxNnodes
write(fileUnit,'(A16,I10)') 'Nnodes', mesh_Nnodes
write(fileUnit,'(A16,I10)') 'NcpElems', mesh_NcpElems
do e = 1_pInt,mesh_NcpElems
t = mesh_element(2,e)
write(fileUnit,'(20(I10))') mesh_element(1_pInt:4_pInt+FE_Nnodes(t),e)
enddo
write(fileUnit,'(A16,I10)') 'Ncellnodes', mesh_Ncellnodes
do n = 1_pInt,mesh_Ncellnodes
write(fileUnit,'(2(I10))') mesh_cellnodeParent(1:2,n)
enddo
write(fileUnit,'(A16,I10)') 'Ncells', mesh_Ncells
do e = 1_pInt,mesh_NcpElems
t = mesh_element(2,e)
g = FE_geomtype(t)
c = FE_celltype(g)
do i = 1_pInt,FE_Nips(g)
write(fileUnit,'(8(I10))') &
mesh_cell(1_pInt:FE_NcellnodesPerCell(c),i,e)
enddo
enddo
close(fileUnit)
end subroutine mesh_write_meshfile
!--------------------------------------------------------------------------------------------------
!> @brief reads mesh description file
!--------------------------------------------------------------------------------------------------
integer function mesh_read_meshfile(filepath)
implicit none
character(len=*), intent(in) :: filepath
integer(pInt), parameter :: fileUnit = 223_pInt
integer(pInt) :: e,i,t,g,n
open(fileUnit,status='old',err=100,iostat=mesh_read_meshfile,action='read',file=filepath)
read(fileUnit,'(TR16,E10.3)',err=100,iostat=mesh_read_meshfile) mesh_unitlength
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_maxNcellnodes
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_maxNips
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_maxNnodes
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_Nnodes
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_NcpElems
if (.not. allocated(mesh_element)) allocate(mesh_element(4_pInt+mesh_maxNnodes,mesh_NcpElems))
mesh_element = 0_pInt
do e = 1_pInt,mesh_NcpElems
read(fileUnit,'(20(I10))',err=100,iostat=mesh_read_meshfile) &
mesh_element(:,e)
enddo
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_Ncellnodes
if (.not. allocated(mesh_cellnodeParent)) allocate(mesh_cellnodeParent(2_pInt,mesh_Ncellnodes))
do n = 1_pInt,mesh_Ncellnodes
read(fileUnit,'(2(I10))',err=100,iostat=mesh_read_meshfile) mesh_cellnodeParent(1:2,n)
enddo
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_Ncells
if (.not. allocated(mesh_cell)) allocate(mesh_cell(FE_maxNcellnodesPerCell,mesh_maxNips,mesh_NcpElems))
do e = 1_pInt,mesh_NcpElems
t = mesh_element(2,e)
g = FE_geomtype(t)
do i = 1_pInt,FE_Nips(g)
read(fileUnit,'(8(I10))',err=100,iostat=mesh_read_meshfile) mesh_cell(:,i,e)
enddo
enddo
close(fileUnit)
mesh_read_meshfile = 0 ! successfully read data
100 continue
end function mesh_read_meshfile
!--------------------------------------------------------------------------------------------------
!> @brief initializes mesh data for use in post processing
!--------------------------------------------------------------------------------------------------
integer function mesh_init_postprocessing(filepath)
implicit none
character(len=*), intent(in) :: filepath
call mesh_build_FEdata
mesh_init_postprocessing = mesh_read_meshfile(filepath)
end function mesh_init_postprocessing
!--------------------------------------------------------------------------------------------------
!> @brief returns global variable mesh_Ncellnodes
!--------------------------------------------------------------------------------------------------

202
code/numerics.f90
View File

@ -236,11 +236,9 @@ subroutine numerics_init
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,ierr);CHKERRQ(ierr)
#endif
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- numerics init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(/,a)') ' <<<+- numerics init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
!$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS...
!$ if(gotDAMASK_NUM_THREADS /= 0) then ! could not get number of threads, set it to 1
@ -489,14 +487,8 @@ subroutine numerics_init
close(FILEUNIT)
else fileExists
#ifdef FEM
if (worldrank == 0) then
#endif
write(6,'(a,/)') ' using standard values'
flush(6)
#ifdef FEM
endif
#endif
endif fileExists
#ifdef Spectral
@ -519,128 +511,126 @@ subroutine numerics_init
!--------------------------------------------------------------------------------------------------
! writing parameters to output
mainProcess3: if (worldrank == 0) then
write(6,'(a24,1x,es8.1)') ' relevantStrain: ',relevantStrain
write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance
write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness
write(6,'(a24,1x,i8)') ' iJacoLpresiduum: ',iJacoLpresiduum
write(6,'(a24,1x,es8.1)') ' pert_Fg: ',pert_Fg
write(6,'(a24,1x,i8)') ' pert_method: ',pert_method
write(6,'(a24,1x,i8)') ' nCryst: ',nCryst
write(6,'(a24,1x,es8.1)') ' subStepMinCryst: ',subStepMinCryst
write(6,'(a24,1x,es8.1)') ' subStepSizeCryst: ',subStepSizeCryst
write(6,'(a24,1x,es8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst
write(6,'(a24,1x,i8)') ' nState: ',nState
write(6,'(a24,1x,i8)') ' nStress: ',nStress
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator
write(6,'(a24,1x,L8)') ' timeSyncing: ',numerics_timeSyncing
write(6,'(a24,1x,L8)') ' analytic Jacobian: ',analyticJaco
write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong
write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength
write(6,'(a24,1x,es8.1)') ' relevantStrain: ',relevantStrain
write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance
write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness
write(6,'(a24,1x,i8)') ' iJacoLpresiduum: ',iJacoLpresiduum
write(6,'(a24,1x,es8.1)') ' pert_Fg: ',pert_Fg
write(6,'(a24,1x,i8)') ' pert_method: ',pert_method
write(6,'(a24,1x,i8)') ' nCryst: ',nCryst
write(6,'(a24,1x,es8.1)') ' subStepMinCryst: ',subStepMinCryst
write(6,'(a24,1x,es8.1)') ' subStepSizeCryst: ',subStepSizeCryst
write(6,'(a24,1x,es8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst
write(6,'(a24,1x,i8)') ' nState: ',nState
write(6,'(a24,1x,i8)') ' nStress: ',nStress
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator
write(6,'(a24,1x,L8)') ' timeSyncing: ',numerics_timeSyncing
write(6,'(a24,1x,L8)') ' analytic Jacobian: ',analyticJaco
write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong
write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength
write(6,'(a24,1x,i8)') ' nHomog: ',nHomog
write(6,'(a24,1x,es8.1)') ' subStepMinHomog: ',subStepMinHomog
write(6,'(a24,1x,es8.1)') ' subStepSizeHomog: ',subStepSizeHomog
write(6,'(a24,1x,es8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog
write(6,'(a24,1x,i8,/)') ' nMPstate: ',nMPstate
write(6,'(a24,1x,i8)') ' nHomog: ',nHomog
write(6,'(a24,1x,es8.1)') ' subStepMinHomog: ',subStepMinHomog
write(6,'(a24,1x,es8.1)') ' subStepSizeHomog: ',subStepSizeHomog
write(6,'(a24,1x,es8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog
write(6,'(a24,1x,i8,/)') ' nMPstate: ',nMPstate
!--------------------------------------------------------------------------------------------------
! RGC parameters
write(6,'(a24,1x,es8.1)') ' aTol_RGC: ',absTol_RGC
write(6,'(a24,1x,es8.1)') ' rTol_RGC: ',relTol_RGC
write(6,'(a24,1x,es8.1)') ' aMax_RGC: ',absMax_RGC
write(6,'(a24,1x,es8.1)') ' rMax_RGC: ',relMax_RGC
write(6,'(a24,1x,es8.1)') ' perturbPenalty_RGC: ',pPert_RGC
write(6,'(a24,1x,es8.1)') ' relevantMismatch_RGC: ',xSmoo_RGC
write(6,'(a24,1x,es8.1)') ' viscosityrate_RGC: ',viscPower_RGC
write(6,'(a24,1x,es8.1)') ' viscositymodulus_RGC: ',viscModus_RGC
write(6,'(a24,1x,es8.1)') ' maxrelaxation_RGC: ',maxdRelax_RGC
write(6,'(a24,1x,es8.1)') ' maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC
write(6,'(a24,1x,es8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC
write(6,'(a24,1x,es8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC
write(6,'(a24,1x,es8.1)') ' aTol_RGC: ',absTol_RGC
write(6,'(a24,1x,es8.1)') ' rTol_RGC: ',relTol_RGC
write(6,'(a24,1x,es8.1)') ' aMax_RGC: ',absMax_RGC
write(6,'(a24,1x,es8.1)') ' rMax_RGC: ',relMax_RGC
write(6,'(a24,1x,es8.1)') ' perturbPenalty_RGC: ',pPert_RGC
write(6,'(a24,1x,es8.1)') ' relevantMismatch_RGC: ',xSmoo_RGC
write(6,'(a24,1x,es8.1)') ' viscosityrate_RGC: ',viscPower_RGC
write(6,'(a24,1x,es8.1)') ' viscositymodulus_RGC: ',viscModus_RGC
write(6,'(a24,1x,es8.1)') ' maxrelaxation_RGC: ',maxdRelax_RGC
write(6,'(a24,1x,es8.1)') ' maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC
write(6,'(a24,1x,es8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC
write(6,'(a24,1x,es8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC
!--------------------------------------------------------------------------------------------------
! Random seeding parameter
write(6,'(a24,1x,i16,/)') ' fixed_seed: ',fixedSeed
if (fixedSeed <= 0_pInt) &
write(6,'(a,/)') ' No fixed Seed: Random is random!'
write(6,'(a24,1x,i16,/)') ' fixed_seed: ',fixedSeed
if (fixedSeed <= 0_pInt) &
write(6,'(a,/)') ' No fixed Seed: Random is random!'
!--------------------------------------------------------------------------------------------------
! gradient parameter
write(6,'(a24,1x,es8.1)') ' charLength: ',charLength
write(6,'(a24,1x,es8.1)') ' residualStiffness: ',residualStiffness
write(6,'(a24,1x,es8.1)') ' charLength: ',charLength
write(6,'(a24,1x,es8.1)') ' residualStiffness: ',residualStiffness
!--------------------------------------------------------------------------------------------------
! openMP parameter
!$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt
!$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt
!--------------------------------------------------------------------------------------------------
! field parameters
write(6,'(a24,1x,i8)') ' itmax: ',itmax
write(6,'(a24,1x,i8)') ' itmin: ',itmin
write(6,'(a24,1x,i8)') ' maxCutBack: ',maxCutBack
write(6,'(a24,1x,i8)') ' maxStaggeredIter: ',stagItMax
write(6,'(a24,1x,i8)') ' vacancyPolyOrder: ',vacancyPolyOrder
write(6,'(a24,1x,i8)') ' hydrogenPolyOrder: ',hydrogenPolyOrder
write(6,'(a24,1x,es8.1)') ' err_struct_tolAbs: ',err_struct_tolAbs
write(6,'(a24,1x,es8.1)') ' err_struct_tolRel: ',err_struct_tolRel
write(6,'(a24,1x,es8.1)') ' err_thermal_tolabs: ',err_thermal_tolabs
write(6,'(a24,1x,es8.1)') ' err_thermal_tolrel: ',err_thermal_tolrel
write(6,'(a24,1x,es8.1)') ' err_damage_tolabs: ',err_damage_tolabs
write(6,'(a24,1x,es8.1)') ' err_damage_tolrel: ',err_damage_tolrel
write(6,'(a24,1x,es8.1)') ' err_vacancyflux_tolabs: ',err_vacancyflux_tolabs
write(6,'(a24,1x,es8.1)') ' err_vacancyflux_tolrel: ',err_vacancyflux_tolrel
write(6,'(a24,1x,es8.1)') ' err_porosity_tolabs: ',err_porosity_tolabs
write(6,'(a24,1x,es8.1)') ' err_porosity_tolrel: ',err_porosity_tolrel
write(6,'(a24,1x,es8.1)') ' err_hydrogenflux_tolabs:',err_hydrogenflux_tolabs
write(6,'(a24,1x,es8.1)') ' err_hydrogenflux_tolrel:',err_hydrogenflux_tolrel
write(6,'(a24,1x,es8.1)') ' vacancyBoundPenalty: ',vacancyBoundPenalty
write(6,'(a24,1x,es8.1)') ' hydrogenBoundPenalty: ',hydrogenBoundPenalty
write(6,'(a24,1x,i8)') ' itmax: ',itmax
write(6,'(a24,1x,i8)') ' itmin: ',itmin
write(6,'(a24,1x,i8)') ' maxCutBack: ',maxCutBack
write(6,'(a24,1x,i8)') ' maxStaggeredIter: ',stagItMax
write(6,'(a24,1x,i8)') ' vacancyPolyOrder: ',vacancyPolyOrder
write(6,'(a24,1x,i8)') ' hydrogenPolyOrder: ',hydrogenPolyOrder
write(6,'(a24,1x,es8.1)') ' err_struct_tolAbs: ',err_struct_tolAbs
write(6,'(a24,1x,es8.1)') ' err_struct_tolRel: ',err_struct_tolRel
write(6,'(a24,1x,es8.1)') ' err_thermal_tolabs: ',err_thermal_tolabs
write(6,'(a24,1x,es8.1)') ' err_thermal_tolrel: ',err_thermal_tolrel
write(6,'(a24,1x,es8.1)') ' err_damage_tolabs: ',err_damage_tolabs
write(6,'(a24,1x,es8.1)') ' err_damage_tolrel: ',err_damage_tolrel
write(6,'(a24,1x,es8.1)') ' err_vacancyflux_tolabs: ',err_vacancyflux_tolabs
write(6,'(a24,1x,es8.1)') ' err_vacancyflux_tolrel: ',err_vacancyflux_tolrel
write(6,'(a24,1x,es8.1)') ' err_porosity_tolabs: ',err_porosity_tolabs
write(6,'(a24,1x,es8.1)') ' err_porosity_tolrel: ',err_porosity_tolrel
write(6,'(a24,1x,es8.1)') ' err_hydrogenflux_tolabs:',err_hydrogenflux_tolabs
write(6,'(a24,1x,es8.1)') ' err_hydrogenflux_tolrel:',err_hydrogenflux_tolrel
write(6,'(a24,1x,es8.1)') ' vacancyBoundPenalty: ',vacancyBoundPenalty
write(6,'(a24,1x,es8.1)') ' hydrogenBoundPenalty: ',hydrogenBoundPenalty
!--------------------------------------------------------------------------------------------------
! spectral parameters
#ifdef Spectral
write(6,'(a24,1x,i8)') ' continueCalculation: ',continueCalculation
write(6,'(a24,1x,L8)') ' memory_efficient: ',memory_efficient
write(6,'(a24,1x,i8)') ' divergence_correction: ',divergence_correction
write(6,'(a24,1x,a)') ' spectral_derivative: ',trim(spectral_derivative)
if(fftw_timelimit<0.0_pReal) then
write(6,'(a24,1x,L8)') ' fftw_timelimit: ',.false.
else
write(6,'(a24,1x,es8.1)') ' fftw_timelimit: ',fftw_timelimit
endif
write(6,'(a24,1x,a)') ' fftw_plan_mode: ',trim(fftw_plan_mode)
write(6,'(a24,1x,i8)') ' fftw_planner_flag: ',fftw_planner_flag
write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma
write(6,'(a24,1x,es8.1)') ' err_stress_tolAbs: ',err_stress_tolAbs
write(6,'(a24,1x,es8.1)') ' err_stress_tolRel: ',err_stress_tolRel
write(6,'(a24,1x,es8.1)') ' err_div_tolAbs: ',err_div_tolAbs
write(6,'(a24,1x,es8.1)') ' err_div_tolRel: ',err_div_tolRel
write(6,'(a24,1x,es8.1)') ' err_curl_tolAbs: ',err_curl_tolAbs
write(6,'(a24,1x,es8.1)') ' err_curl_tolRel: ',err_curl_tolRel
write(6,'(a24,1x,es8.1)') ' polarAlpha: ',polarAlpha
write(6,'(a24,1x,es8.1)') ' polarBeta: ',polarBeta
write(6,'(a24,1x,a)') ' spectral solver: ',trim(spectral_solver)
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_defaultOptions)//' '//trim(petsc_options)
write(6,'(a24,1x,i8)') ' continueCalculation: ',continueCalculation
write(6,'(a24,1x,L8)') ' memory_efficient: ',memory_efficient
write(6,'(a24,1x,i8)') ' divergence_correction: ',divergence_correction
write(6,'(a24,1x,a)') ' spectral_derivative: ',trim(spectral_derivative)
if(fftw_timelimit<0.0_pReal) then
write(6,'(a24,1x,L8)') ' fftw_timelimit: ',.false.
else
write(6,'(a24,1x,es8.1)') ' fftw_timelimit: ',fftw_timelimit
endif
write(6,'(a24,1x,a)') ' fftw_plan_mode: ',trim(fftw_plan_mode)
write(6,'(a24,1x,i8)') ' fftw_planner_flag: ',fftw_planner_flag
write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma
write(6,'(a24,1x,es8.1)') ' err_stress_tolAbs: ',err_stress_tolAbs
write(6,'(a24,1x,es8.1)') ' err_stress_tolRel: ',err_stress_tolRel
write(6,'(a24,1x,es8.1)') ' err_div_tolAbs: ',err_div_tolAbs
write(6,'(a24,1x,es8.1)') ' err_div_tolRel: ',err_div_tolRel
write(6,'(a24,1x,es8.1)') ' err_curl_tolAbs: ',err_curl_tolAbs
write(6,'(a24,1x,es8.1)') ' err_curl_tolRel: ',err_curl_tolRel
write(6,'(a24,1x,es8.1)') ' polarAlpha: ',polarAlpha
write(6,'(a24,1x,es8.1)') ' polarBeta: ',polarBeta
write(6,'(a24,1x,a)') ' spectral solver: ',trim(spectral_solver)
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_defaultOptions)//' '//trim(petsc_options)
#endif
!--------------------------------------------------------------------------------------------------
! spectral parameters
#ifdef FEM
write(6,'(a24,1x,i8)') ' integrationOrder: ',integrationOrder
write(6,'(a24,1x,i8)') ' structOrder: ',structOrder
write(6,'(a24,1x,i8)') ' thermalOrder: ',thermalOrder
write(6,'(a24,1x,i8)') ' damageOrder: ',damageOrder
write(6,'(a24,1x,i8)') ' vacancyfluxOrder: ',vacancyfluxOrder
write(6,'(a24,1x,i8)') ' porosityOrder: ',porosityOrder
write(6,'(a24,1x,i8)') ' hydrogenfluxOrder: ',hydrogenfluxOrder
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_defaultOptions)//' '//trim(petsc_options)
write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation
write(6,'(a24,1x,i8)') ' integrationOrder: ',integrationOrder
write(6,'(a24,1x,i8)') ' structOrder: ',structOrder
write(6,'(a24,1x,i8)') ' thermalOrder: ',thermalOrder
write(6,'(a24,1x,i8)') ' damageOrder: ',damageOrder
write(6,'(a24,1x,i8)') ' vacancyfluxOrder: ',vacancyfluxOrder
write(6,'(a24,1x,i8)') ' porosityOrder: ',porosityOrder
write(6,'(a24,1x,i8)') ' hydrogenfluxOrder: ',hydrogenfluxOrder
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_defaultOptions)//' '//trim(petsc_options)
write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation
#endif
endif mainProcess3
!--------------------------------------------------------------------------------------------------

29
code/prec.f90
View File

@ -130,30 +130,17 @@ subroutine prec_init
iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
implicit none
integer(pInt) :: worldrank = 0_pInt
#ifdef PETSc
#include <petsc/finclude/petscsys.h>
PetscErrorCode :: ierr
#endif
external :: &
quit, &
MPI_Comm_rank, &
MPI_Abort
#ifdef PETSc
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
#endif
quit
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- prec init -+>>>'
write(6,'(/,a)') ' <<<+- prec init -+>>>'
#include "compilation_info.f90"
write(6,'(a,i3)') ' Bytes for pReal: ',pReal
write(6,'(a,i3)') ' Bytes for pInt: ',pInt
write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt
write(6,'(a,e10.3)') ' NaN: ', DAMASK_NaN
write(6,'(a,l3)') ' NaN != NaN: ',DAMASK_NaN /= DAMASK_NaN
write(6,'(a,l3,/)') ' NaN check passed ',prec_isNAN(DAMASK_NaN)
endif mainProcess
write(6,'(a,i3)') ' Bytes for pReal: ',pReal
write(6,'(a,i3)') ' Bytes for pInt: ',pInt
write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt
write(6,'(a,e10.3)') ' NaN: ', DAMASK_NaN
write(6,'(a,l3)') ' NaN != NaN: ',DAMASK_NaN /= DAMASK_NaN
write(6,'(a,l3,/)') ' NaN check passed ',prec_isNAN(DAMASK_NaN)
if ((.not. prec_isNaN(DAMASK_NaN)) .or. (DAMASK_NaN == DAMASK_NaN)) call quit(9000)
realloc_lhs_test = [1_pInt,2_pInt]

55
code/spectral_damage.f90
View File

@ -42,7 +42,6 @@ module spectral_damage
integer(pInt), private :: totalIter = 0_pInt !< total iteration in current increment
real(pReal), dimension(3,3), private :: D_ref
real(pReal), private :: mobility_ref
character(len=1024), private :: incInfo
public :: &
spectral_damage_init, &
@ -50,21 +49,7 @@ module spectral_damage
spectral_damage_forward, &
spectral_damage_destroy
external :: &
VecDestroy, &
DMDestroy, &
DMDACreate3D, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
PETScFinalize, &
SNESDestroy, &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber, &
SNESSolve, &
SNESSetDM, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions, &
SNESCreate, &
MPI_Abort, &
MPI_Bcast, &
MPI_Allreduce
@ -90,15 +75,30 @@ subroutine spectral_damage_init()
damage_nonlocal_getMobility
implicit none
integer(pInt), dimension(:), allocatable :: localK
integer(pInt) :: proc
integer(pInt) :: i, j, k, cell
DM :: damage_grid
Vec :: uBound, lBound
PetscErrorCode :: ierr
PetscObject :: dummy
integer(pInt), dimension(:), allocatable :: localK
integer(pInt) :: proc
integer(pInt) :: i, j, k, cell
character(len=100) :: snes_type
external :: &
SNESCreate, &
SNESSetOptionsPrefix, &
DMDACreate3D, &
SNESSetDM, &
DMDAGetCorners, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
SNESSetFromOptions, &
SNESGetType, &
VecSet, &
DMGetGlobalVector, &
DMRestoreGlobalVector, &
SNESVISetVariableBounds
mainProcess: if (worldrank == 0_pInt) then
write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -194,12 +194,18 @@ type(tSolutionState) function spectral_damage_solution(guess,timeinc,timeinc_old
integer(pInt) :: i, j, k, cell
PetscInt ::position
PetscReal :: minDamage, maxDamage, stagNorm, solnNorm
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscErrorCode :: ierr
SNESConvergedReason :: reason
external :: &
VecMin, &
VecMax, &
SNESSolve, &
SNESGetConvergedReason
spectral_damage_solution%converged =.false.
!--------------------------------------------------------------------------------------------------
@ -353,10 +359,13 @@ subroutine spectral_damage_forward(guess,timeinc,timeinc_old,loadCaseTime)
timeinc, &
loadCaseTime !< remaining time of current load case
logical, intent(in) :: guess
PetscErrorCode :: ierr
integer(pInt) :: i, j, k, cell
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
external :: &
SNESGetDM
if (cutBack) then
damage_current = damage_lastInc
@ -400,6 +409,10 @@ subroutine spectral_damage_destroy()
implicit none
PetscErrorCode :: ierr
external :: &
VecDestroy, &
SNESDestroy
call VecDestroy(solution,ierr); CHKERRQ(ierr)
call SNESDestroy(damage_snes,ierr); CHKERRQ(ierr)

193
code/spectral_interface.f90
View File

@ -13,9 +13,7 @@ module DAMASK_interface
pInt
implicit none
private
#ifdef PETSc
#include <petsc/finclude/petscsys.h>
#endif
logical, public, protected :: appendToOutFile = .false. !< Append to existing spectralOut file (in case of restart, not in case of regridding)
integer(pInt), public, protected :: spectralRestartInc = 1_pInt !< Increment at which calculation starts
character(len=1024), public, protected :: &
@ -44,13 +42,10 @@ contains
!> @brief initializes the solver by interpreting the command line arguments. Also writes
!! information on computation to screen
!--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
subroutine DAMASK_interface_init()
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
implicit none
character(len=1024), optional, intent(in) :: &
loadCaseParameterIn, & !< if using the f2py variant, the -l argument of DAMASK_spectral.exe
geometryParameterIn !< if using the f2py variant, the -g argument of DAMASK_spectral.exe
character(len=1024) :: &
commandLine, & !< command line call as string
loadCaseArg ='', & !< -l argument given to DAMASK_spectral.exe
@ -67,9 +62,7 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
chunkPos
integer, dimension(8) :: &
dateAndTime ! type default integer
#ifdef PETSc
PetscErrorCode :: ierr
#endif
external :: &
quit,&
MPI_Comm_rank,&
@ -77,9 +70,10 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
MPI_Init_Thread, &
MPI_abort
open(6, encoding='UTF-8') ! for special characters in output
!--------------------------------------------------------------------------------------------------
! PETSc Init
#ifdef PETSc
#ifdef _OPENMP
call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,ierr);CHKERRQ(ierr) ! in case of OpenMP, don't rely on PETScInitialize doing MPI init
if (threadLevel<MPI_THREAD_FUNNELED) then
@ -89,104 +83,95 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
#endif
call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code
CHKERRQ(ierr) ! this is a macro definition, it is case sensitive
open(6, encoding='UTF-8')
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
#endif
mainProcess: if (worldrank == 0) then
call date_and_time(values = dateAndTime)
write(6,'(/,a)') ' <<<+- DAMASK_spectral -+>>>'
write(6,'(/,a)') ' Version: '//DAMASKVERSION
write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
#include "compilation_info.f90"
if (output_unit /= 6) then
write(output_unit,'(a)') 'STDOUT != 6'
call quit(1_pInt)
endif
else mainProcess
close(6) ! disable output for non-master processes (open 6 to rank specific file for debug)
open(6,file='/dev/null',status='replace') ! close(6) alone will leave some temp files in cwd
endif mainProcess
if ( present(loadcaseParameterIn) .and. present(geometryParameterIn)) then ! both mandatory parameters given in function call
geometryArg = geometryParameterIn
loadcaseArg = loadcaseParameterIn
commandLine = 'n/a'
else if ( .not.( present(loadcaseParameterIn) .and. present(geometryParameterIn))) then ! none parameters given in function call, trying to get them from command line
call get_command(commandLine)
chunkPos = IIO_stringPos(commandLine)
do i = 1, chunkPos(1)
tag = IIO_lc(IIO_stringValue(commandLine,chunkPos,i)) ! extract key
select case(tag)
case ('-h','--help')
mainProcess2: if (worldrank == 0) then
write(6,'(a)') ' #######################################################################'
write(6,'(a)') ' DAMASK_spectral:'
write(6,'(a)') ' The spectral method boundary value problem solver for'
write(6,'(a)') ' the Düsseldorf Advanced Material Simulation Kit'
write(6,'(a,/)')' #######################################################################'
write(6,'(a,/)')' Valid command line switches:'
write(6,'(a)') ' --geom (-g, --geometry)'
write(6,'(a)') ' --load (-l, --loadcase)'
write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory, -d, --directory)'
write(6,'(a)') ' --restart (-r, --rs)'
write(6,'(a)') ' --regrid (--rg)'
write(6,'(a)') ' --help (-h)'
write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Mandatory arguments:'
write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom.geom'
write(6,'(a)') ' Specifies the location of the geometry definition file,'
write(6,'(a)') ' if no extension is given, .geom will be appended.'
write(6,'(a)') ' "PathToGeomFile" will be the working directory if not specified'
write(6,'(a)') ' via --workingdir.'
write(6,'(a)') ' Make sure the file "material.config" exists in the working'
write(6,'(a)') ' directory.'
write(6,'(a)') ' For further configuration place "numerics.config"'
write(6,'(a)')' and "numerics.config" in that directory.'
write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile.load'
write(6,'(a)') ' Specifies the location of the load case definition file,'
write(6,'(a)') ' if no extension is given, .load will be appended.'
write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Optional arguments:'
write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory'
write(6,'(a)') ' Specifies the working directory and overwrites the default'
write(6,'(a)') ' "PathToGeomFile".'
write(6,'(a)') ' Make sure the file "material.config" exists in the working'
write(6,'(a)') ' directory.'
write(6,'(a)') ' For further configuration place "numerics.config"'
write(6,'(a)')' and "numerics.config" in that directory.'
write(6,'(/,a)')' --restart XX'
write(6,'(a)') ' Reads in total increment No. XX-1 and continues to'
write(6,'(a)') ' calculate total increment No. XX.'
write(6,'(a)') ' Appends to existing results file '
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".'
write(6,'(a)') ' Works only if the restart information for total increment'
write(6,'(a)') ' No. XX-1 is available in the working directory.'
write(6,'(/,a)')' --regrid XX'
write(6,'(a)') ' Reads in total increment No. XX-1 and continues to'
write(6,'(a)') ' calculate total increment No. XX.'
write(6,'(a)') ' Attention: Overwrites existing results file '
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".'
write(6,'(a)') ' Works only if the restart information for total increment'
write(6,'(a)') ' No. XX-1 is available in the working directory.'
write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Help:'
write(6,'(/,a)')' --help'
write(6,'(a,/)')' Prints this message and exits'
call quit(0_pInt) ! normal Termination
endif mainProcess2
case ('-l', '--load', '--loadcase')
loadcaseArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt)
case ('-g', '--geom', '--geometry')
geometryArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt)
case ('-w', '-d', '--wd', '--directory', '--workingdir', '--workingdirectory')
workingDirArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt)
case ('-r', '--rs', '--restart')
spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt)
appendToOutFile = .true.
case ('--rg', '--regrid')
spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt)
appendToOutFile = .false.
end select
enddo
endif
call date_and_time(values = dateAndTime)
write(6,'(/,a)') ' <<<+- DAMASK_spectral -+>>>'
write(6,'(/,a)') ' Version: '//DAMASKVERSION
write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',&
dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
#include "compilation_info.f90"
call get_command(commandLine)
chunkPos = IIO_stringPos(commandLine)
do i = 1, chunkPos(1)
tag = IIO_lc(IIO_stringValue(commandLine,chunkPos,i)) ! extract key
select case(tag)
case ('-h','--help')
mainProcess2: if (worldrank == 0) then
write(6,'(a)') ' #######################################################################'
write(6,'(a)') ' DAMASK_spectral:'
write(6,'(a)') ' The spectral method boundary value problem solver for'
write(6,'(a)') ' the Düsseldorf Advanced Material Simulation Kit'
write(6,'(a,/)')' #######################################################################'
write(6,'(a,/)')' Valid command line switches:'
write(6,'(a)') ' --geom (-g, --geometry)'
write(6,'(a)') ' --load (-l, --loadcase)'
write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory, -d, --directory)'
write(6,'(a)') ' --restart (-r, --rs)'
write(6,'(a)') ' --help (-h)'
write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Mandatory arguments:'
write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom.geom'
write(6,'(a)') ' Specifies the location of the geometry definition file,'
write(6,'(a)') ' if no extension is given, .geom will be appended.'
write(6,'(a)') ' "PathToGeomFile" will be the working directory if not specified'
write(6,'(a)') ' via --workingdir.'
write(6,'(a)') ' Make sure the file "material.config" exists in the working'
write(6,'(a)') ' directory.'
write(6,'(a)') ' For further configuration place "numerics.config"'
write(6,'(a)')' and "numerics.config" in that directory.'
write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile.load'
write(6,'(a)') ' Specifies the location of the load case definition file,'
write(6,'(a)') ' if no extension is given, .load will be appended.'
write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Optional arguments:'
write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory'
write(6,'(a)') ' Specifies the working directory and overwrites the default'
write(6,'(a)') ' "PathToGeomFile".'
write(6,'(a)') ' Make sure the file "material.config" exists in the working'
write(6,'(a)') ' directory.'
write(6,'(a)') ' For further configuration place "numerics.config"'
write(6,'(a)')' and "numerics.config" in that directory.'
write(6,'(/,a)')' --restart XX'
write(6,'(a)') ' Reads in total increment No. XX-1 and continues to'
write(6,'(a)') ' calculate total increment No. XX.'
write(6,'(a)') ' Appends to existing results file '
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".'
write(6,'(a)') ' Works only if the restart information for total increment'
write(6,'(a)') ' No. XX-1 is available in the working directory.'
write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Help:'
write(6,'(/,a)')' --help'
write(6,'(a,/)')' Prints this message and exits'
call quit(0_pInt) ! normal Termination
endif mainProcess2
case ('-l', '--load', '--loadcase')
loadcaseArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt)
case ('-g', '--geom', '--geometry')
geometryArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt)
case ('-w', '-d', '--wd', '--directory', '--workingdir', '--workingdirectory')
workingDirArg = IIO_stringValue(commandLine,chunkPos,i+1_pInt)
case ('-r', '--rs', '--restart')
spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt)
appendToOutFile = .true.
end select
enddo
if (len(trim(loadcaseArg)) == 0 .or. len(trim(geometryArg)) == 0) then
write(6,'(a)') ' Please specify geometry AND load case (-h for help)'

67
code/spectral_mech_AL.f90
View File

@ -22,7 +22,7 @@ module spectral_mech_AL
DAMASK_spectral_solverAL_label = 'al'
!--------------------------------------------------------------------------------------------------
! derived types
! derived types
type(tSolutionParams), private :: params
real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
@ -31,7 +31,7 @@ module spectral_mech_AL
DM, private :: da
SNES, private :: snes
Vec, private :: solution_vec
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
@ -72,21 +72,7 @@ module spectral_mech_AL
AL_forward, &
AL_destroy
external :: &
VecDestroy, &
DMDestroy, &
DMDACreate3D, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
PETScFinalize, &
SNESDestroy, &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber, &
SNESSolve, &
SNESSetDM, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions, &
SNESCreate, &
MPI_Abort, &
MPI_Bcast, &
MPI_Allreduce
@ -136,11 +122,22 @@ subroutine AL_init
integer(pInt) :: proc
character(len=1024) :: rankStr
if (worldrank == 0_pInt) then
external :: &
SNESCreate, &
SNESSetOptionsPrefix, &
DMDACreate3D, &
SNESSetDM, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions
mainProcess: if (worldrank == 0_pInt) then
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif
endif mainProcess
!--------------------------------------------------------------------------------------------------
! allocate global fields
@ -150,7 +147,7 @@ subroutine AL_init
allocate (F_lambdaDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! PETSc Init
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3
@ -185,10 +182,10 @@ subroutine AL_init
'reading values of increment ', restartInc - 1_pInt, ' from file'
flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_realFile(777,'F'//trim(rankStr), trim(getSolverJobName()),size(F))
call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F))
read (777,rec=1) F
close (777)
call IO_read_realFile(777,'F_lastInc'//trim(rankStr), trim(getSolverJobName()),size(F_lastInc))
call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc
close (777)
call IO_read_realFile(777,'F_lambda'//trim(rankStr),trim(getSolverJobName()),size(F_lambda))
@ -214,15 +211,14 @@ subroutine AL_init
F_lambda_lastInc = F_lastInc
endif restart
call Utilities_updateIPcoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), &
0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3)
nullify(F)
nullify(F_lambda)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc
readRestart: if (restartInc > 1_pInt) then
restartRead: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file'
@ -236,7 +232,7 @@ subroutine AL_init
call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
read (777,rec=1) C_minMaxAvg
close (777)
endif readRestart
endif restartRead
call Utilities_updateGamma(C_minMaxAvg,.True.)
C_scale = C_minMaxAvg
@ -263,7 +259,7 @@ type(tSolutionState) function &
use FEsolving, only: &
restartWrite, &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
@ -286,6 +282,10 @@ type(tSolutionState) function &
PetscErrorCode :: ierr
SNESConvergedReason :: reason
external :: &
SNESSolve, &
SNESGetConvergedReason
incInfo = incInfoIn
!--------------------------------------------------------------------------------------------------
@ -298,7 +298,7 @@ type(tSolutionState) function &
endif
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
! set module wide availabe data
mask_stress = P_BC%maskFloat
params%P_BC = P_BC%values
params%rotation_BC = rotation_BC
@ -387,6 +387,10 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
integer(pInt) :: &
i, j, k, e
external :: &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber
F => x_scal(1:3,1:3,1,&
XG_RANGE,YG_RANGE,ZG_RANGE)
F_lambda => x_scal(1:3,1:3,2,&
@ -414,7 +418,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', &
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
math_transpose33(F_aim)
math_transpose33(F_aim)
flush(6)
endif
endif newIteration
@ -507,7 +511,7 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr
fnorm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode ::ierr
PetscErrorCode :: ierr
real(pReal) :: &
curlTol, &
divTol, &
@ -704,6 +708,11 @@ subroutine AL_destroy()
implicit none
PetscErrorCode :: ierr
external :: &
VecDestroy, &
SNESDestroy, &
DMDestroy
call VecDestroy(solution_vec,ierr); CHKERRQ(ierr)
call SNESDestroy(snes,ierr); CHKERRQ(ierr)
call DMDestroy(da,ierr); CHKERRQ(ierr)

95
code/spectral_mech_Basic.f90
View File

@ -48,7 +48,7 @@ module spectral_mech_basic
C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness
S = 0.0_pReal !< current compliance (filled up with zeros)
S = 0.0_pReal !< current compliance (filled up with zeros)
real(pReal), private :: err_stress, err_div
logical, private :: ForwardData
integer(pInt), private :: &
@ -61,21 +61,7 @@ module spectral_mech_basic
BasicPETSc_forward, &
basicPETSc_destroy
external :: &
VecDestroy, &
DMDestroy, &
DMDACreate3D, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
PETScFinalize, &
SNESDestroy, &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber, &
SNESSolve, &
SNESSetDM, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions, &
SNESCreate, &
MPI_Abort, &
MPI_Bcast, &
MPI_Allreduce
@ -105,7 +91,7 @@ subroutine basicPETSc_init
use spectral_utilities, only: &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
utilities_updateIPcoords, &
Utilities_updateIPcoords, &
wgt
use mesh, only: &
grid, &
@ -115,15 +101,28 @@ subroutine basicPETSc_init
implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
PetscScalar, dimension(:,:,:,:), pointer :: F
PetscErrorCode :: ierr
PetscObject :: dummy
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
PetscErrorCode :: ierr
PetscObject :: dummy
PetscScalar, pointer, dimension(:,:,:,:) :: F
integer(pInt), dimension(:), allocatable :: localK
integer(pInt) :: proc
character(len=1024) :: rankStr
external :: &
SNESCreate, &
SNESSetOptionsPrefix, &
DMDACreate3D, &
SNESSetDM, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions
mainProcess: if (worldrank == 0_pInt) then
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -147,9 +146,9 @@ subroutine basicPETSc_init
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & ! global grid
1, 1, worldsize, &
1 , 1, worldsize, &
9, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
grid (1),grid (2),localK, & ! local grid
grid(1),grid(2),localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
@ -195,10 +194,9 @@ subroutine basicPETSc_init
temp33_Real, &
.false., &
math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc
restartRead: if (restartInc > 1_pInt) then
restartRead: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file'
@ -243,19 +241,24 @@ type(tSolutionState) function &
timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment
loadCaseTime !< remaining time of current load case
logical, intent(in) :: &
guess
type(tBoundaryCondition), intent(in) :: &
P_BC, &
F_BC
character(len=*), intent(in) :: &
incInfoIn
real(pReal), dimension(3,3), intent(in) :: rotation_BC
logical, intent(in) :: &
guess
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscErrorCode :: ierr
SNESConvergedReason :: reason
external :: &
SNESSolve, &
SNESGetConvergedReason
incInfo = incInfoIn
!--------------------------------------------------------------------------------------------------
@ -263,9 +266,9 @@ type(tSolutionState) function &
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg)
if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite)
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
! set module wide availabe data
mask_stress = P_BC%maskFloat
params%P_BC = P_BC%values
params%rotation_BC = rotation_BC
@ -292,7 +295,7 @@ end function BasicPETSc_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the AL residual vector
!> @brief forms the basic residual vector
!--------------------------------------------------------------------------------------------------
subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
use numerics, only: &
@ -312,10 +315,11 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
debug_spectral, &
debug_spectralRotation
use spectral_utilities, only: &
wgt, &
tensorField_real, &
utilities_FFTtensorForward, &
utilities_FFTtensorBackward, &
utilities_fourierGammaConvolution, &
utilities_FFTtensorBackward, &
Utilities_constitutiveResponse, &
Utilities_divergenceRMS
use IO, only: &
@ -338,11 +342,15 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
PetscObject :: dummy
PetscErrorCode :: ierr
external :: &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
newIteration: if (totalIter <= PETScIter) then
newIteration: if(totalIter <= PETScIter) then
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
totalIter = totalIter + 1_pInt
@ -351,7 +359,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', &
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
math_transpose33(F_aim)
flush(6)
@ -401,7 +409,7 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
worldrank
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
@ -415,10 +423,10 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
real(pReal) :: &
divTol, &
stressTol
divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs)
stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)
converged: if ((totalIter >= itmin .and. &
all([ err_div/divTol, &
err_stress/stressTol ] < 1.0_pReal)) &
@ -451,21 +459,21 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use numerics, only: &
worldrank
use mesh, only: &
grid, &
grid3
use spectral_utilities, only: &
Utilities_calculateRate, &
Utilities_forwardField, &
utilities_updateIPcoords, &
Utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use IO, only: &
IO_write_JobRealFile
use FEsolving, only: &
restartWrite
use numerics, only: &
worldrank
implicit none
real(pReal), intent(in) :: &
@ -478,8 +486,9 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
real(pReal), dimension(3,3), intent(in) :: rotation_BC
logical, intent(in) :: &
guess
PetscErrorCode :: ierr
PetscScalar, pointer :: F(:,:,:,:)
PetscErrorCode :: ierr
character(len=1024) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
@ -508,7 +517,7 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
write (777,rec=1) C_volAvgLastInc
close(777)
endif
endif
endif
call utilities_updateIPcoords(F)
@ -538,6 +547,7 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
endif
F_aim = F_aim + f_aimDot * timeinc
!--------------------------------------------------------------------------------------------------
@ -558,6 +568,11 @@ subroutine BasicPETSc_destroy()
implicit none
PetscErrorCode :: ierr
external :: &
VecDestroy, &
SNESDestroy, &
DMDestroy
call VecDestroy(solution_vec,ierr); CHKERRQ(ierr)
call SNESDestroy(snes,ierr); CHKERRQ(ierr)
call DMDestroy(da,ierr); CHKERRQ(ierr)

67
code/spectral_mech_Polarisation.f90
View File

@ -22,7 +22,7 @@ module spectral_mech_Polarisation
DAMASK_spectral_solverPolarisation_label = 'polarisation'
!--------------------------------------------------------------------------------------------------
! derived types
! derived types
type(tSolutionParams), private :: params
real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
@ -31,7 +31,7 @@ module spectral_mech_Polarisation
DM, private :: da
SNES, private :: snes
Vec, private :: solution_vec
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
@ -57,7 +57,7 @@ module spectral_mech_Polarisation
S = 0.0_pReal, & !< current compliance (filled up with zeros)
C_scale = 0.0_pReal, &
S_scale = 0.0_pReal
real(pReal), private :: &
err_BC, & !< deviation from stress BC
err_curl, & !< RMS of curl of F
@ -72,21 +72,7 @@ module spectral_mech_Polarisation
Polarisation_forward, &
Polarisation_destroy
external :: &
VecDestroy, &
DMDestroy, &
DMDACreate3D, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
PETScFinalize, &
SNESDestroy, &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber, &
SNESSolve, &
SNESSetDM, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions, &
SNESCreate, &
MPI_Abort, &
MPI_Bcast, &
MPI_Allreduce
@ -136,11 +122,22 @@ subroutine Polarisation_init
integer(pInt) :: proc
character(len=1024) :: rankStr
if (worldrank == 0_pInt) then
external :: &
SNESCreate, &
SNESSetOptionsPrefix, &
DMDACreate3D, &
SNESSetDM, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions
mainProcess: if (worldrank == 0_pInt) then
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif
endif mainProcess
!--------------------------------------------------------------------------------------------------
! allocate global fields
@ -150,7 +147,7 @@ subroutine Polarisation_init
allocate (F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! PETSc Init
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3
@ -163,7 +160,7 @@ subroutine Polarisation_init
grid(1),grid(2),grid(3), & ! global grid
1 , 1, worldsize, &
18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
grid (1),grid (2),localK, & ! local grid
grid(1),grid(2),localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
@ -182,7 +179,7 @@ subroutine Polarisation_init
restart: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment', restartInc - 1_pInt, 'from file'
'reading values of increment ', restartInc - 1_pInt, ' from file'
flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F))
@ -221,7 +218,7 @@ subroutine Polarisation_init
nullify(F_tau)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc
readRestart: if (restartInc > 1_pInt) then
restartRead: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file'
@ -235,7 +232,7 @@ subroutine Polarisation_init
call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
read (777,rec=1) C_minMaxAvg
close (777)
endif readRestart
endif restartRead
call Utilities_updateGamma(C_minMaxAvg,.True.)
C_scale = C_minMaxAvg
@ -262,7 +259,7 @@ type(tSolutionState) function &
use FEsolving, only: &
restartWrite, &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
@ -285,6 +282,10 @@ type(tSolutionState) function &
PetscErrorCode :: ierr
SNESConvergedReason :: reason
external :: &
SNESSolve, &
SNESGetConvergedReason
incInfo = incInfoIn
!--------------------------------------------------------------------------------------------------
@ -385,7 +386,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
PetscErrorCode :: ierr
integer(pInt) :: &
i, j, k, e
external :: &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber
F => x_scal(1:3,1:3,1,&
XG_RANGE,YG_RANGE,ZG_RANGE)
F_tau => x_scal(1:3,1:3,2,&
@ -505,7 +510,7 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
fnorm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode ::ierr
PetscErrorCode :: ierr
real(pReal) :: &
curlTol, &
divTol, &
@ -631,7 +636,8 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
write (777,rec=1) C_volAvgLastInc
close(777)
endif
endif
endif
call utilities_updateIPcoords(F)
if (cutBack) then
@ -701,6 +707,11 @@ subroutine Polarisation_destroy()
implicit none
PetscErrorCode :: ierr
external :: &
VecDestroy, &
SNESDestroy, &
DMDestroy
call VecDestroy(solution_vec,ierr); CHKERRQ(ierr)
call SNESDestroy(snes,ierr); CHKERRQ(ierr)
call DMDestroy(da,ierr); CHKERRQ(ierr)

50
code/spectral_thermal.f90
View File

@ -42,7 +42,6 @@ module spectral_thermal
integer(pInt), private :: totalIter = 0_pInt !< total iteration in current increment
real(pReal), dimension(3,3), private :: D_ref
real(pReal), private :: mobility_ref
character(len=1024), private :: incInfo
public :: &
spectral_thermal_init, &
@ -50,21 +49,7 @@ module spectral_thermal
spectral_thermal_forward, &
spectral_thermal_destroy
external :: &
VecDestroy, &
DMDestroy, &
DMDACreate3D, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
PETScFinalize, &
SNESDestroy, &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber, &
SNESSolve, &
SNESSetDM, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions, &
SNESCreate, &
MPI_Abort, &
MPI_Bcast, &
MPI_Allreduce
@ -99,10 +84,20 @@ subroutine spectral_thermal_init
integer(pInt) :: proc
integer(pInt) :: i, j, k, cell
DM :: thermal_grid
PetscScalar, pointer :: x_scal(:,:,:)
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
PetscObject :: dummy
external :: &
SNESCreate, &
SNESSetOptionsPrefix, &
DMDACreate3D, &
SNESSetDM, &
DMDAGetCorners, &
DMCreateGlobalVector, &
DMDASNESSetFunctionLocal, &
SNESSetFromOptions
mainProcess: if (worldrank == 0_pInt) then
write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -154,6 +149,8 @@ subroutine spectral_thermal_init
x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current
call DMDAVecRestoreArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! thermal reference diffusion update
cell = 0_pInt
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
@ -171,7 +168,7 @@ subroutine spectral_thermal_init
end subroutine spectral_thermal_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the Basic PETSC scheme with internal iterations
!> @brief solution for the spectral thermal scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_old,loadCaseTime)
use numerics, only: &
@ -196,12 +193,18 @@ type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_ol
integer(pInt) :: i, j, k, cell
PetscInt :: position
PetscReal :: minTemperature, maxTemperature, stagNorm, solnNorm
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscErrorCode :: ierr
SNESConvergedReason :: reason
external :: &
VecMin, &
VecMax, &
SNESSolve, &
SNESGetConvergedReason
spectral_thermal_solution%converged =.false.
!--------------------------------------------------------------------------------------------------
@ -355,8 +358,11 @@ subroutine spectral_thermal_forward(guess,timeinc,timeinc_old,loadCaseTime)
logical, intent(in) :: guess
integer(pInt) :: i, j, k, cell
DM :: dm_local
PetscScalar, pointer :: x_scal(:,:,:)
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
external :: &
SNESGetDM
if (cutBack) then
temperature_current = temperature_lastInc
@ -405,6 +411,10 @@ subroutine spectral_thermal_destroy()
implicit none
PetscErrorCode :: ierr
external :: &
VecDestroy, &
SNESDestroy
call VecDestroy(solution,ierr); CHKERRQ(ierr)
call SNESDestroy(thermal_snes,ierr); CHKERRQ(ierr)

159
configure vendored
View File

@ -22,19 +22,6 @@ class extendableOption(Option):
Option.take_action(self, action, dest, opt, value, values, parser)
# -----------------------------
def filePresent(paths,files,warning=False):
for path in paths:
for file in files:
if os.path.isfile(os.path.join(path,file)): return True
if warning: print "Warning: %s not found in %s"%(','.join(files),','.join(paths))
return False
########################################################
# MAIN
########################################################
@ -42,35 +29,15 @@ def filePresent(paths,files,warning=False):
parser = OptionParser(option_class=extendableOption, usage='%prog options', description = """
Configures the compilation and installation of DAMASK
""" + string.replace('$Id$','\n','\\n')
)
#--- determine default compiler ----------------------------------------------------------------------
compiler = os.getenv('F90')
if compiler == None:
compiler = 'ifort' if subprocess.call(['which', 'ifort'], stdout=subprocess.PIPE, stderr=subprocess.PIPE) == 0 \
else 'gfortran'
#--- default option values --------------------------------------------------------------------------
BLAS_order = ['IMKL','ACML','LAPACK','OPENBLAS']
""")
defaults={'DAMASK_BIN':'depending on access rights',
'F90':compiler,
'FFTW_ROOT':'/usr',
'MSC_ROOT' :'/msc',
'DAMASK_NUM_THREADS':4,
'MARC_VERSION':'2015',
'blasType':'LAPACK',
'blasRoot':{'LAPACK' :'/usr',
'ACML' :'/opt/acml6.1.0',
'IMKL' : os.getenv('MKLROOT') if os.getenv('MKLROOT') else '/opt/intel/composerxe/mkl',
'OPENBLAS' :'/usr',
},
'spectralOptions':{},
}
#--- if local config file exists, read, otherwise assume global config file ------------------------
configFile = os.path.join(os.getenv('HOME'),'.damask/damask.conf') \
if os.path.isfile(os.path.join(os.getenv('HOME'),'.damask/damask.conf')) \
@ -91,129 +58,25 @@ try:
defaults['DAMASK_NUM_THREADS'] = int(value)
if key == 'DAMASK_BIN':
defaults['DAMASK_BIN'] = value
if key in ['F90','FFTW_ROOT','MSC_ROOT','spectralOptions','MARC_VERSION']:
defaults[key] = value
for theKey in reversed(BLAS_order):
if key == theKey+'_ROOT' and value != None and value != '':
defaults['blasType'] = theKey
defaults['blasRoot'][theKey] = value
except IOError:
pass
parser.add_option('--prefix', dest='prefix', metavar='string',
help='location of (links to) DAMASK executables [%default]')
parser.add_option('--with-FC','--with-fc',
dest='compiler', metavar='string',
help='F90 compiler [%default]')
parser.add_option('--with-FFTW-dir','--with-fftw-dir',
dest='fftwRoot', metavar='string',
help='root directory of FFTW [%default]')
parser.add_option('--with-MSC-dir','--with-msc-dir',
dest='mscRoot', metavar='string',
help='root directory of MSC.Marc/Mentat [%default]')
parser.add_option('--with-MARC-version','--with-marc-version',
dest='marcVersion', metavar='string',
help='version of MSC.Marc/Mentat [%default]')
parser.add_option('--with-OMP-threads','--with-omp-threads',
dest='threads', type='int', metavar='int',
help='number of openMP threads [%default]')
parser.add_option('--with-BLAS-type','--with-blas-type',
dest='blasType', metavar='string',
help='type of BLAS/LAPACK library [%default] {{{}}}'.format(','.join(BLAS_order)))
parser.add_option('--with-BLAS-dir','--with-blas-dir',
dest='blasRoot',metavar='string',
help='root directory of BLAS/LAPACK library [%default]')
parser.add_option('--with-spectral-options', dest='spectraloptions', action='extend', metavar='<string LIST>',
help='options for compilation of spectral solver')
parser.set_defaults(prefix = defaults['DAMASK_BIN'])
parser.set_defaults(compiler = defaults['F90'])
parser.set_defaults(fftwRoot = defaults['FFTW_ROOT'])
parser.set_defaults(mscRoot = defaults['MSC_ROOT'])
parser.set_defaults(marcVersion = defaults['MARC_VERSION'])
parser.set_defaults(threads = defaults['DAMASK_NUM_THREADS'])
parser.set_defaults(blasType = defaults['blasType'])
#--- set default for blasRoot depending on current option (or default) for blasType --------------------
blasType = defaults['blasType'].upper()
for i, arg in enumerate(sys.argv):
if arg.lower().startswith('--with-blas-type'):
if arg.lower().endswith('--with-blas-type'):
blasType = sys.argv[i+1].upper()
else:
blasType = sys.argv[i][17:].upper()
if blasType not in BLAS_order:
blasType = defaults['blasType'].upper()
parser.set_defaults(blasRoot = defaults['blasRoot'][blasType])
parser.set_defaults(spectraloptions = [])
(options,filenames) = parser.parse_args()
#--- consistency checks --------------------------------------------------------------------------------
options.compiler = options.compiler.lower()
options.blasType = options.blasType.upper()
options.fftwRoot = os.path.normpath(options.fftwRoot)
options.mscRoot = os.path.normpath(options.mscRoot)
options.blasRoot = os.path.normpath(options.blasRoot)
locations = {
'FFTW' : [os.path.join(options.fftwRoot,'lib64'),os.path.join(options.fftwRoot,'lib')],
'LAPACK' : [os.path.join(options.blasRoot,'lib64'),os.path.join(options.blasRoot,'lib')],
'OPENBLAS': [os.path.join(options.blasRoot,'lib64'),os.path.join(options.blasRoot,'lib')],
'ACML' : [os.path.join(options.blasRoot,'%s64/lib'%options.compiler)],
'ACML_mp' : [os.path.join(options.blasRoot,'%s64_mp/lib'%options.compiler)],
'IMKL' : [os.path.join(options.blasRoot,'lib/intel64')],
}
libraries = {
'FFTW' : [
'libfftw3.so','libfftw3.a',
'libfftw3_threads.so','libfftw3_threads.a',
],
'LAPACK' : [
'liblapack.so','liblapack.a','liblapack.dylib',
],
'OPENBLAS' : [
'libopenblas.so','libopenblas.a','libopenblas.dylib',
],
'ACML' : [
'libacml.so','libacml.a',
],
'ACML_mp' : [
'libacml_mp.so','libacml_mp.a',
],
'IMKL' : [
'libmkl_core.so','libmkl_core.a',
'libmkl_sequential.so','libmkl_sequential.a',
'libmkl_intel_thread.so','libmkl_intel_thread.a',
'libmkl_intel_lp64.so','libmkl_intel_lp64.a',
'libmkl_gnu_thread.so','libmkl_gnu_thread.a',
'libmkl_gf_lp64.so','libmkl_gf_lp64.a',
],
}
if options.compiler not in ['ifort','gfortran']:
print('Error: Unknown compiler option: %s'%options.compiler)
sys.exit(1)
if not subprocess.call(['which', options.compiler], stdout=subprocess.PIPE, stderr=subprocess.PIPE) == 0:
print('Compiler Warning: executable %s not found!'%options.compiler)
if not os.path.isdir(options.mscRoot):
print('Warning: MSC root directory %s not found!'%options.mscRoot)
filePresent(locations['FFTW'],libraries['FFTW'],warning=True)
if options.blasType in ['LAPACK','OPENBLAS','IMKL']:
filePresent(locations[options.blasType],libraries[options.blasType],warning=True)
elif options.blasType == 'ACML':
filePresent(locations[options.blasType],libraries[options.blasType],warning=True)
filePresent(locations[options.blasType+'_mp'],libraries[options.blasType+'_mp'],warning=True)
else:
print('Error: Unknown BLAS/LAPACK library: %s'%options.blasType)
sys.exit(1)
#--- read config file if present to keep comments and order ---------------------------------------
output = []
@ -228,12 +91,6 @@ try:
if items[0] == 'DAMASK_BIN':
line = '%s=%s'%(items[0],options.prefix)
options.prefix ='depending on access rights'
if items[0] == 'F90':
line = '%s=%s'%(items[0],options.compiler)
options.compiler =''
if items[0] == 'FFTW_ROOT':
line = '%s=%s'%(items[0],options.fftwRoot)
options.fftwRoot =''
if items[0] == 'MSC_ROOT':
line = '%s=%s'%(items[0],options.mscRoot)
options.mscRoot =''
@ -243,14 +100,6 @@ try:
if items[0] == 'DAMASK_NUM_THREADS':
line = '%s=%s'%(items[0],options.threads)
options.threads =''
for blasType in defaults['blasRoot'].keys():
if items[0] == '%s_ROOT'%blasType and items[0] == '%s_ROOT'%options.blasType:
line = '%s=%s'%(items[0],options.blasRoot)
options.blasType=''
elif items[0] == '#%s_ROOT'%blasType and items[0] == '#%s_ROOT'%options.blasType:
line = '%s=%s'%(items[0][1:],options.blasRoot)
options.blasType=''
elif items[0] == '%s_ROOT'%blasType: line = '#'+line
for spectralOption in options.spectraloptions:
[key,value] = re.split('[= ]',spectralOption)[0:2]
if key == items[0]:
@ -264,18 +113,12 @@ except IOError:
for opt, value in options.__dict__.items():
if opt == 'prefix' and value != 'depending on access rights':
output.append('DAMASK_BIN=%s'%value)
if opt == 'compiler' and value != '':
output.append('F90=%s'%value)
if opt == 'fftwRoot' and value != '':
output.append('FFTW_ROOT=%s'%value)
if opt == 'mscRoot' and value != '':
output.append('MSC_ROOT=%s'%value)
if opt == 'marcVersion' and value != '':
output.append('MARC_VERSION=%s'%value)
if opt == 'threads' and value != '':
output.append('DAMASK_NUM_THREADS=%s'%value)
if opt == 'blasType' and value != '':
output.append('%s_ROOT=%s'%(options.blasType,options.blasRoot))
for spectralOption in options.spectraloptions:
output.append(spectralOption)

141
installation/compile_CoreModule.py
View File

@ -1,141 +0,0 @@
#!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*-
import os,sys,glob,subprocess,shlex
from damask import Environment
from damask import version as DAMASKVERSION
# compiles fortran code for Python
scriptID = '$Id$'
damaskEnv = Environment()
baseDir = damaskEnv.relPath('installation/')
codeDir = damaskEnv.relPath('code/')
keywords=['IMKL_ROOT','ACML_ROOT','LAPACK_ROOT','FFTW_ROOT','F90']
options={}
#--- getting options from damask.conf or, if not present, from envinronment -----------------------
for option in keywords:
try:
value = damaskEnv.options[option]
except:
value = os.getenv(option)
if value is None: value = '' # env not set
options[option]=value
#--- overwrite default options with keyword=value pair from argument list to mimic make behavior --
for i, arg in enumerate(sys.argv):
for option in keywords:
if arg.startswith(option):
options[option] = sys.argv[i][len(option)+1:]
#--- check for valid compiler and set options -----------------------------------------------------
compilers = ['ifort','gfortran']
if options['F90'] not in compilers:
sys.exit('compiler "F90" (in installation/options or as Shell variable) has to be one out of: %s'%(', '.join(compilers)))
compiler = {
'gfortran': '--fcompiler=gnu95 --f90flags="-fPIC -fno-range-check -xf95-cpp-input -std=f2008 -fall-intrinsics'+\
' -fdefault-real-8 -fdefault-double-8"',
'ifort': '--fcompiler=intelem --f90flags="-fPIC -fpp -stand f08 -diag-disable 5268 -assume byterecl'+\
' -real-size 64 -integer-size 32 -shared-intel"',
}[options['F90']]
#--- option not depending on compiler -------------------------------------------------------------
compileOptions = ' -DSpectral -DFLOAT=8 -DINT=4 -I%s/lib -DDAMASKVERSION=\\\\\"\\\"%s\\\\\"\\\"'%(damaskEnv.rootDir(),DAMASKVERSION)
#--- this saves the path of libraries to core.so, hence it is known during runtime ----------------
if options['F90'] == 'gfortran':
# solved error: Undefined symbols for architecture x86_64: "_PyArg_ParseTupleAndKeywords"
# as found on https://lists.macosforge.org/pipermail/macports-dev/2013-May/022735.html
LDFLAGS = '-shared -Wl,-undefined,dynamic_lookup'
else:
# some f2py versions/configurations compile with openMP, so linking against openMP is needed
# to prevent errors during loading of core module
LDFLAGS = ' -openmp -Wl'
#--- run path of for fftw during runtime ----------------------------------------------------------
LDFLAGS += ',-rpath,%s/lib,-rpath,%s/lib64'%(options['FFTW_ROOT'],options['FFTW_ROOT'])
# see http://cens.ioc.ee/pipermail/f2py-users/2003-December/000621.html
if options['IMKL_ROOT']:
if options['F90'] == 'gfortran':
arch = 'gf'
elif options['F90'] == 'ifort':
arch = 'intel'
lib_lapack = '-L%s/lib/intel64 -lmkl_%s_lp64 -lmkl_core -lmkl_sequential -lm'\
%(options['IMKL_ROOT'],arch)
LDFLAGS +=',-rpath,%s/lib/intel64'%(options['IMKL_ROOT'])
elif options['ACML_ROOT'] != '':
lib_lapack = '-L%s/%s64/lib -lacml'%(options['ACML_ROOT'],options['F90'])
LDFLAGS +=',-rpath,%s/%s64/lib'%(options['ACML_ROOT'],options['F90'])
elif options['LAPACK_ROOT'] != '':
lib_lapack = '-L%s/lib -L%s/lib64 -llapack'%(options['LAPACK_ROOT'],options['LAPACK_ROOT'])
LDFLAGS +=',-rpath,%s/lib,-rpath,%s/lib64'%(options['LAPACK_ROOT'],options['LAPACK_ROOT'])
#--------------------------------------------------------------------------------------------------
# f2py does not (yet) support setting of special flags for the linker, hence they must be set via
# environment variable ----------------------------------------------------------------------------
my_env = os.environ
my_env["LDFLAGS"] = LDFLAGS
#--- delete old file ------------------------------------------------------------------------------
try:
os.remove(os.path.join(damaskEnv.relPath('lib/damask'),'core.so'))
except OSError, e:
print ("Error when deleting: %s - %s." % (e.filename,e.strerror))
# The following command is used to compile the fortran files and make the functions defined
# in damask.core.pyf available for python in the module core.so
# It uses the fortran wrapper f2py that is included in the numpy package to construct the
# module core.so out of the fortran code in the f90 files
# For the generation of the pyf file use the following lines:
###########################################################################
#'f2py -h damask.core.pyf' +\
#' --overwrite-signature --no-lower prec.f90 DAMASK_spectral_interface.f90 math.f90 mesh.f90,...'
###########################################################################
os.chdir(codeDir)
cmd = 'f2py damask.core.pyf' +\
' -c --no-lower %s'%(compiler) +\
compileOptions+\
' C_routines.c'+\
' system_routines.f90'+\
' prec.f90'+\
' spectral_interface.f90'+\
' IO.f90'+\
' libs.f90'+\
' numerics.f90'+\
' debug.f90'+\
' math.f90'+\
' FEsolving.f90'+\
' mesh.f90'+\
' core_quit.f90'+\
' -L%s/lib -lfftw3'%(options['FFTW_ROOT'])+\
' %s'%lib_lapack
print('Executing: '+cmd)
try:
subprocess.call(shlex.split(cmd),env=my_env)
except subprocess.CalledProcessError:
print('build failed')
except OSError:
print ('f2py not found')
try:
os.rename(os.path.join(codeDir,'core.so'),\
os.path.join(damaskEnv.relPath('lib/damask'),'core.so'))
except:
pass
modules = glob.glob('*.mod')
for module in modules:
print 'removing', module
os.remove(module)
#--- check if compilation of core module was successful -------------------------------------------
try:
with open(damaskEnv.relPath('lib/damask/core.so')) as f: pass
except IOError as e:
print '*********\n* core.so not found, compilation of core modules was not successful\n*********'

1230
lib/IR_Precision.f90
View File

File diff suppressed because it is too large Load Diff

909
lib/Lib_Base64.f90
View File

@ -1,909 +0,0 @@
!> @ingroup Library
!> @{
!> @defgroup Lib_Base64Library Lib_Base64
!> base64 encoding/decoding library
!> @}
!> @ingroup Interface
!> @{
!> @defgroup Lib_Base64Interface Lib_Base64
!> base64 encoding/decoding library
!> @}
!> @ingroup PublicProcedure
!> @{
!> @defgroup Lib_Base64PublicProcedure Lib_Base64
!> base64 encoding/decoding library
!> @}
!> @ingroup PrivateProcedure
!> @{
!> @defgroup Lib_Base64PrivateProcedure Lib_Base64
!> base64 encoding/decoding library
!> @}
!> @ingroup GlobalVarPar
!> @{
!> @defgroup Lib_Base64GlobalVarPar Lib_Base64
!> base64 encoding/decoding library
!> @}
!> @ingroup PrivateVarPar
!> @{
!> @defgroup Lib_Base64PrivateVarPar Lib_Base64
!> base64 encoding/decoding library
!> @}
!> This module contains base64 encoding/decoding procedures.
!> @todo \b Decoding: Implement decoding functions.
!> @todo \b DocComplete: Complete the documentation.
!> @ingroup Lib_Base64Library
module Lib_Base64
!-----------------------------------------------------------------------------------------------------------------------------------
USE IR_Precision ! Integers and reals precision definition.
!-----------------------------------------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------------------------------------
implicit none
private
public:: b64_encode
!public:: b64_decode
public:: pack_data
!-----------------------------------------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------------------------------------
!> @ingroup Lib_Base64GlobalVarPar
!> @{
!> @}
!> @ingroup Lib_Base64PrivateVarPar
!> @{
character(64):: base64="ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/" !< Base64 alphabet.
!> @}
!-----------------------------------------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------------------------------------
!> @brief Subroutine for encoding numbers (integer and real) to base64.
!> @ingroup Lib_Base64Interface
interface b64_encode
module procedure b64_encode_R8_a, &
b64_encode_R4_a, &
b64_encode_I8_a, &
b64_encode_I4_a, &
b64_encode_I2_a, &
b64_encode_I1_a
endinterface
!!> @brief Subroutine for decoding numbers (integer and real) from base64.
!!> @ingroup Lib_Base64Interface
!interface b64_decode
! module procedure b64_decode_R8_a, &
! b64_decode_R4_a, &
! b64_decode_I8_a, &
! b64_decode_I4_a, &
! b64_decode_I2_a, &
! b64_decode_I1_a
!endinterface
!> @brief Subroutine for packing different kinds of data into single I1P array. This is useful for encoding different kinds
!> variables into a single stream of bits.
!> @ingroup Lib_Base64Interface
interface pack_data
module procedure pack_data_R8_R4,pack_data_R8_I8,pack_data_R8_I4,pack_data_R8_I2,pack_data_R8_I1, &
pack_data_R4_R8,pack_data_R4_I8,pack_data_R4_I4,pack_data_R4_I2,pack_data_R4_I1, &
pack_data_I8_R8,pack_data_I8_R4,pack_data_I8_I4,pack_data_I8_I2,pack_data_I8_I1, &
pack_data_I4_R8,pack_data_I4_R4,pack_data_I4_I8,pack_data_I4_I2,pack_data_I4_I1, &
pack_data_I2_R8,pack_data_I2_R4,pack_data_I2_I8,pack_data_I2_I4,pack_data_I2_I1, &
pack_data_I1_R8,pack_data_I1_R4,pack_data_I1_I8,pack_data_I1_I4,pack_data_I1_I2
endinterface
!-----------------------------------------------------------------------------------------------------------------------------------
contains
!> @ingroup Lib_Base64PrivateProcedure
!> @{
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R8_R4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R8P), intent(IN):: a1(1:) !< Firs data stream.
real(R4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R8_R4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R8_I8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R8P), intent(IN):: a1(1:) !< First data stream.
integer(I8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R8_I8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R8_I4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R8P), intent(IN):: a1(1:) !< First data stream.
integer(I4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R8_I4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R8_I2(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R8P), intent(IN):: a1(1:) !< First data stream.
integer(I2P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R8_I2
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R8_I1(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R8P), intent(IN):: a1(1:) !< First data stream.
integer(I1P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R8_I1
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R4_R8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R4P), intent(IN):: a1(1:) !< Firs data stream.
real(R8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R4_R8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R4_I8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R4P), intent(IN):: a1(1:) !< First data stream.
integer(I8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R4_I8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R4_I4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R4P), intent(IN):: a1(1:) !< First data stream.
integer(I4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R4_I4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R4_I2(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R4P), intent(IN):: a1(1:) !< First data stream.
integer(I2P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R4_I2
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R4_I1(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R4P), intent(IN):: a1(1:) !< First data stream.
integer(I1P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R4_I1
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I8_R8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I8P), intent(IN):: a1(1:) !< First data stream.
real(R8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I8_R8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I8_R4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I8P), intent(IN):: a1(1:) !< First data stream.
real(R4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I8_R4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I8_I4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I8P), intent(IN):: a1(1:) !< First data stream.
integer(I4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I8_I4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I8_I2(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I8P), intent(IN):: a1(1:) !< First data stream.
integer(I2P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I8_I2
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I8_I1(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I8P), intent(IN):: a1(1:) !< First data stream.
integer(I1P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I8_I1
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I4_R8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: a1(1:) !< First data stream.
real(R8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I4_R8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I4_R4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: a1(1:) !< First data stream.
real(R4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I4_R4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I4_I8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: a1(1:) !< First data stream.
integer(I8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I4_I8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I4_I2(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: a1(1:) !< First data stream.
integer(I2P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I4_I2
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I4_I1(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: a1(1:) !< First data stream.
integer(I1P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I4_I1
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I2_R8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I2P), intent(IN):: a1(1:) !< First data stream.
real(R8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I2_R8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I2_R4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I2P), intent(IN):: a1(1:) !< First data stream.
real(R4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I2_R4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I2_I8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I2P), intent(IN):: a1(1:) !< First data stream.
integer(I8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I2_I8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I2_I4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I2P), intent(IN):: a1(1:) !< First data stream.
integer(I4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I2_I4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I2_I1(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I2P), intent(IN):: a1(1:) !< First data stream.
integer(I1P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I2_I1
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I1_R8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: a1(1:) !< First data stream.
real(R8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I1_R8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I1_R4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: a1(1:) !< First data stream.
real(R4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I1_R4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I1_I8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: a1(1:) !< First data stream.
integer(I8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I1_I8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I1_I4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: a1(1:) !< First data stream.
integer(I4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I1_I4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I1_I2(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: a1(1:) !< First data stream.
integer(I2P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I1_I2
!> @brief Subroutine for encoding bits (must be multiple of 24 bits) into base64 charcaters code (of length multiple of 4).
!> @note The bits stream are encoded in chunks of 24 bits as the following example (in little endian order):
!> @code
!> +--first octet--+-second octet--+--third octet--+
!> |7 6 5 4 3 2 1 0|7 6 5 4 3 2 1 0|7 6 5 4 3 2 1 0|
!> +-----------+---+-------+-------+---+-----------+
!> |5 4 3 2 1 0|5 4 3 2 1 0|5 4 3 2 1 0|5 4 3 2 1 0|
!> +--1.index--+--2.index--+--3.index--+--4.index--+
!> @endcode
!> The 4 indexes are stored into 4 elements 8 bits array, thus 2 bits of each array element are not used.
pure subroutine encode_bits(bits,padd,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: bits(1:) !< Bits to be encoded.
integer(I4P), intent(IN):: padd !< Number of padding characters ('=').
character(1), intent(OUT):: code(1:) !< Characters code.
integer(I1P):: sixb(1:4) !< 6 bits slices (stored into 8 bits integer) of 24 bits input.
integer(I8P):: c,e !< Counters.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
c = 1_I8P
do e=1_I8P,size(bits,dim=1,kind=I8P),3_I8P ! loop over array elements: 3 bytes (24 bits) scanning
sixb = 0_I1P
call mvbits(bits(e ),2,6,sixb(1),0)
call mvbits(bits(e ),0,2,sixb(2),4) ; call mvbits(bits(e+1),4,4,sixb(2),0)
call mvbits(bits(e+1),0,4,sixb(3),2) ; call mvbits(bits(e+2),6,2,sixb(3),0)
call mvbits(bits(e+2),0,6,sixb(4),0)
sixb = sixb + 1_I1P
code(c :c )(1:1) = base64(sixb(1):sixb(1))
code(c+1:c+1)(1:1) = base64(sixb(2):sixb(2))
code(c+2:c+2)(1:1) = base64(sixb(3):sixb(3))
code(c+3:c+3)(1:1) = base64(sixb(4):sixb(4))
c = c + 4_I8P
enddo
if (padd>0) code(size(code,dim=1)-padd+1:)(1:1)='='
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine encode_bits
!> @brief Subroutine for encoding array numbers to base64 (R8P).
pure subroutine b64_encode_R8_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
real(R8P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_R8_a
!> @brief Subroutine for encoding array numbers to base64 (R4P).
pure subroutine b64_encode_R4_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
real(R4P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_R4_a
!> @brief Subroutine for encoding array numbers to base64 (I8P).
pure subroutine b64_encode_I8_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
integer(I8P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_I8_a
!> @brief Subroutine for encoding array numbers to base64 (I4P).
pure subroutine b64_encode_I4_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
integer(I4P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_I4_a
!> @brief Subroutine for encoding array numbers to base64 (I2P).
pure subroutine b64_encode_I2_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
integer(I2P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_I2_a
!> @brief Subroutine for encoding array numbers to base64 (I1P).
pure subroutine b64_encode_I1_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
integer(I1P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_I1_a
!!> @brief Subroutine for decoding array numbers from base64 (R8P).
!pure subroutine b64_decode_R8_a(code,n)
!!--------------------------------------------------------------------------------------------------------------------------------
!implicit none
!real(R8P), intent(OUT):: n(1:) !< Number to be decoded.
!character(ncR8P*size(n,dim=1)), intent(IN):: code !< Encoded number.
!integer(I4P):: c,d !< Counters.
!!--------------------------------------------------------------------------------------------------------------------------------
!!--------------------------------------------------------------------------------------------------------------------------------
!d = 1_I4P
!do c=1,len(code),ncR8P
! call b64_decode_R8_s(code=code(c:c+ncR8P-1),n=n(d))
! d = d + 1_I4P
!enddo
!return
!!--------------------------------------------------------------------------------------------------------------------------------
!endsubroutine b64_decode_R8_a
!> @}
endmodule Lib_Base64

6070
lib/Lib_VTK_IO.f90
View File

File diff suppressed because it is too large Load Diff

1
lib/damask/.gitignore vendored
View File

@ -1,2 +1,3 @@
core.so
corientation.so
*.pyx

45
lib/damask/__init__.py
View File

@ -1,7 +1,7 @@
# -*- coding: UTF-8 no BOM -*-
"""Main aggregator"""
import sys, os
import os
with open(os.path.join(os.path.dirname(__file__),'../../VERSION')) as f:
version = f.readline()[:-1]
@ -10,47 +10,14 @@ from .environment import Environment # noqa
from .asciitable import ASCIItable # noqa
from .config import Material # noqa
from .colormaps import Colormap, Color # noqa
from .orientation import Quaternion, Rodrigues, Symmetry, Orientation # noqa
# try:
# from .corientation import Quaternion, Rodrigues, Symmetry, Orientation
# print "Import Cython version of Orientation module"
# except:
# from .orientation import Quaternion, Rodrigues, Symmetry, Orientation
try:
from .corientation import Quaternion, Rodrigues, Symmetry, Orientation # noqa
print "Import Cython version of Orientation module"
except:
from .orientation import Quaternion, Rodrigues, Symmetry, Orientation # noqa
#from .block import Block # only one class
from .result import Result # noqa
from .geometry import Geometry # noqa
from .solver import Solver # noqa
from .test import Test # noqa
from .util import extendableOption # noqa
try:
from . import core
# cleaning up namespace
###################################################################################################
# capitalize according to convention
core.IO = core.io
core.FEsolving = core.fesolving
core.DAMASK_interface = core.damask_interface
# remove modulePrefix_
core.prec.init = core.prec.prec_init
core.DAMASK_interface.init = core.DAMASK_interface.DAMASK_interface_init
core.IO.init = core.IO.IO_init
core.numerics.init = core.numerics.numerics_init
core.debug.init = core.debug.debug_init
core.math.init = core.math.math_init
core.math.tensorAvg = core.math.math_tensorAvg
core.FEsolving.init = core.FEsolving.FE_init
core.mesh.init = core.mesh.mesh_init
core.mesh.nodesAroundCentres = core.mesh.mesh_nodesAroundCentres
core.mesh.deformedCoordsFFT = core.mesh.mesh_deformedCoordsFFT
core.mesh.volumeMismatch = core.mesh.mesh_volumeMismatch
core.mesh.shapeMismatch = core.mesh.mesh_shapeMismatch
except (ImportError,AttributeError) as e:
core = None # from http://www.python.org/dev/peps/pep-0008/
if os.path.split(sys.argv[0])[1] not in ('symLink_Processing.py',
'compile_CoreModule.py',
):
sys.stderr.write('\nWARNING: Core module (Fortran code) not available, \n'\
"try to run 'make processing'\n"\
'Error message when importing core.so: %s\n\n'%e)

1
lib/damask/colormaps.py
View File

@ -15,6 +15,7 @@ class Color():
__slots__ = [
'model',
'color',
'__dict__',
]

1277
lib/damask/corientation.pyx
View File

File diff suppressed because it is too large Load Diff

2
lib/damask/solver/abaqus.py
View File

@ -13,7 +13,7 @@ class Abaqus(Solver):
import subprocess
process = subprocess.Popen(['abaqus', 'information=release'],stdout = subprocess.PIPE,stderr = subprocess.PIPE)
self.version = process.stdout.readlines()[1].split()[1]
print self.version
print(self.version)
else:
self.version = version

4
lib/damask/util.py
View File

@ -415,7 +415,7 @@ def leastsqBound(func, x0, args=(), bounds=None, Dfun=None, full_output=0,
try:
cov_x = inv(np.dot(np.transpose(R), R))
except LinAlgError as inverror:
print inverror
print(inverror)
pass
return (x, cov_x) + retval[1:-1] + (mesg, info)
else:
@ -464,4 +464,4 @@ def curve_fit_bound(f, xdata, ydata, p0=None, sigma=None, bounds=None, **kw):
else:
pcov = np.inf
return (popt, pcov, infodict, errmsg, ier) if return_full else (popt, pcov)
return (popt, pcov, infodict, errmsg, ier) if return_full else (popt, pcov)

1226
lib/fftw3.f03
View File

File diff suppressed because it is too large Load Diff

433
processing/post/3Dvisualize.py
View File

@ -1,433 +0,0 @@
#!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*-
import os,sys,re,fnmatch,vtk
import numpy as np
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
def outStdout(cmd,locals):
if cmd[0:3] == '(!)':
exec(cmd[3:])
elif cmd[0:3] == '(?)':
cmd = eval(cmd[3:])
print cmd
else:
print cmd
return
def outFile(cmd,locals):
if cmd[0:3] == '(!)':
exec(cmd[3:])
elif cmd[0:3] == '(?)':
cmd = eval(cmd[3:])
locals['filepointer'].write(cmd+'\n')
else:
locals['filepointer'].write(cmd+'\n')
return
def output(cmds,locals,dest):
for cmd in cmds:
if isinstance(cmd,list):
output(cmd,locals,dest)
else:
{\
'File': outFile,\
'Stdout': outStdout,\
}[dest](str(cmd),locals)
return
def transliterateToFloat(x):
try:
return float(x)
except:
return 0.0
def unravel(item):
if hasattr(item,'__contains__'): return ' '.join(map(unravel,item))
else: return str(item)
# ++++++++++++++++++++++++++++++++++++++++++++++++++++
def vtk_writeASCII_mesh(mesh,data,res,sep):
"""function writes data array defined on a hexahedral mesh (geometry)"""
info = {\
'tensor': {'name':'tensor','len':9},\
'vector': {'name':'vector','len':3},\
'scalar': {'name':'scalar','len':1},\
'double': {'name':'scalar','len':2},\
'triple': {'name':'scalar','len':3},\
'quadruple': {'name':'scalar','len':4},\
}
N1 = (res[0]+1)*(res[1]+1)*(res[2]+1)
N = res[0]*res[1]*res[2]
cmds = [\
'# vtk DataFile Version 3.1',
'powered by %s'%scriptID,
'ASCII',
'DATASET UNSTRUCTURED_GRID',
'POINTS %i double'%N1,
[[['\t'.join(map(str,mesh[:,i,j,k])) for i in range(res[0]+1)] for j in range(res[1]+1)] for k in range(res[2]+1)],
'CELLS %i %i'%(N,N*9),
]
# cells
for z in range (res[2]):
for y in range (res[1]):
for x in range (res[0]):
base = z*(res[1]+1)*(res[0]+1)+y*(res[0]+1)+x
cmds.append('8 '+'\t'.join(map(str,[ \
base,
base+1,
base+res[0]+2,
base+res[0]+1,
base+(res[1]+1)*(res[0]+1),
base+(res[1]+1)*(res[0]+1)+1,
base+(res[1]+1)*(res[0]+1)+res[0]+2,
base+(res[1]+1)*(res[0]+1)+res[0]+1,
])))
cmds += [\
'CELL_TYPES %i'%N,
['12']*N,
'CELL_DATA %i'%N,
]
for type in data:
plural = {True:'',False:'S'}[type.lower().endswith('s')]
for item in data[type]['_order_']:
cmds += [\
'%s %s double'%(info[type]['name'].upper()+plural,item),
{True:'LOOKUP_TABLE default',False:''}[info[type]['name'][:3]=='sca'],
[[[sep.join(map(unravel,data[type][item][:,j,k]))] for j in range(res[1])] for k in range(res[2])],
]
return cmds
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def vtk_writeASCII_points(coordinates,data,res,sep):
"""function writes data array defined on a point field"""
N = res[0]*res[1]*res[2]
cmds = [\
'# vtk DataFile Version 3.1',
'powered by %s'%scriptID,
'ASCII',
'DATASET UNSTRUCTURED_GRID',
'POINTS %i double'%N,
[[['\t'.join(map(str,coordinates[i,j,k])) for i in range(res[0])] for j in range(res[1])] for k in range(res[2])],
'CELLS %i %i'%(N,N*2),
['1\t%i'%i for i in range(N)],
'CELL_TYPES %i'%N,
['1']*N,
'POINT_DATA %i'%N,
]
for type in data:
plural = {True:'',False:'S'}[type.lower().endswith('s')]
for item in data[type]:
cmds += [\
'%s %s double'%(type.upper()+plural,item),
{True:'LOOKUP_TABLE default',False:''}[type.lower()[:3]=='sca'],
[[[sep.join(map(unravel,data[type][item][:,j,k]))] for j in range(res[1])] for k in range(res[2])],
]
return cmds
# ----------------------- MAIN -------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog [options] datafile[s]', description = """
Produce VTK file from data field.
Coordinates are taken from (consecutive) x, y, and z columns.
""", version = scriptID)
sepChoices = ['n','t','s']
parser.add_option('-s', '--scalar', dest='scalar', action='extend', metavar = '<string LIST>',
help='list of single scalars to visualize')
parser.add_option( '--double', dest='double', action='extend', metavar = '<string LIST>',
help='list of two scalars to visualize')
parser.add_option( '--triple', dest='triple', action='extend', metavar = '<string LIST>',
help='list of three scalars to visualize')
parser.add_option( '--quadruple', dest='quadruple', action='extend', metavar = '<string LIST>',
help='list of four scalars to visualize')
parser.add_option('-v', '--vector', dest='vector', action='extend', metavar = '<string LIST>',
help='list of vectors to visualize')
parser.add_option('-t', '--tensor', dest='tensor', action='extend', metavar = '<string LIST>',
help='list of tensors to visualize')
parser.add_option('-d', '--deformation', dest='defgrad', metavar = 'string',
help='heading of deformation gradient columns [%default]')
parser.add_option('--reference', dest='undeformed', action='store_true',
help='map results to reference (undeformed) configuration [%default]')
parser.add_option('-c','--cell', dest='cell', action='store_true',
help='data is cell-centered [%default]')
parser.add_option('-p','--vertex', dest='cell', action='store_false',
help='data is vertex-centered')
parser.add_option('--mesh', dest='output_mesh', action='store_true',
help='produce VTK mesh file [%default]')
parser.add_option('--nomesh', dest='output_mesh', action='store_false',
help='omit VTK mesh file')
parser.add_option('--points', dest='output_points', action='store_true',
help='produce VTK points file [%default]')
parser.add_option('--nopoints', dest='output_points', action='store_false',
help='omit VTK points file')
parser.add_option('--separator', dest='separator', type='choice', choices=sepChoices, metavar='string',
help='data separator {%s} [t]'%(' '.join(map(str,sepChoices))))
parser.add_option('--scaling', dest='scaling', action='extend', metavar = '<float LIST>',
help='scaling of fluctuation')
parser.add_option('-u', '--unitlength', dest='unitlength', type='float', metavar = 'float',
help='set unit length for 2D model [%default]')
parser.set_defaults(defgrad = 'f')
parser.set_defaults(separator = 't')
parser.set_defaults(scalar = [])
parser.set_defaults(double = [])
parser.set_defaults(triple = [])
parser.set_defaults(quadruple = [])
parser.set_defaults(vector = [])
parser.set_defaults(tensor = [])
parser.set_defaults(output_mesh = True)
parser.set_defaults(output_points = False)
parser.set_defaults(scaling = [])
parser.set_defaults(undeformed = False)
parser.set_defaults(unitlength = 0.0)
parser.set_defaults(cell = True)
sep = {'n': '\n', 't': '\t', 's': ' '}
(options, args) = parser.parse_args()
options.scaling += [1.0 for i in xrange(max(0,3-len(options.scaling)))]
options.scaling = map(float, options.scaling)
for filename in args:
if not os.path.exists(filename):
continue
file = open(filename)
content = file.readlines()
file.close()
m = re.search('(\d+)\s*head', content[0].lower())
if m is None:
continue
print filename,'\n'
sys.stdout.flush()
headrow = int(m.group(1))
headings = content[headrow].split()
column = {}
matches = {}
maxcol = 0
locol = -1
for col,head in enumerate(headings):
if head == {True:'1_pos',False:'1_nodeinitialcoord'}[options.cell]:
locol = col
maxcol = max(maxcol,col+3)
break
if locol < 0:
print 'missing coordinates..!'
continue
column['tensor'] = {}
matches['tensor'] = {}
for label in [options.defgrad] + options.tensor:
column['tensor'][label] = -1
for col,head in enumerate(headings):
if head == label or head == '1_'+label:
column['tensor'][label] = col
maxcol = max(maxcol,col+9)
matches['tensor'][label] = [label]
break
if not options.undeformed and column['tensor'][options.defgrad] < 0:
print 'missing deformation gradient "%s"..!'%options.defgrad
continue
column['vector'] = {}
matches['vector'] = {}
for label in options.vector:
column['vector'][label] = -1
for col,head in enumerate(headings):
if head == label or head == '1_'+label:
column['vector'][label] = col
maxcol = max(maxcol,col+3)
matches['vector'][label] = [label]
break
for length,what in enumerate(['scalar','double','triple','quadruple']):
column[what] = {}
labels = eval("options.%s"%what)
matches[what] = {}
for col,head in enumerate(headings):
for needle in labels:
if fnmatch.fnmatch(head,needle):
column[what][head] = col
maxcol = max(maxcol,col+1+length)
if needle not in matches[what]:
matches[what][needle] = [head]
else:
matches[what][needle] += [head]
values = np.array(sorted([map(transliterateToFloat,line.split()[:maxcol]) for line in content[headrow+1:]],
key=lambda x:(x[locol+0],x[locol+1],x[locol+2])),'d') # sort with z as fastest and x as slowest index
values2 = np.array([map(transliterateToFloat,line.split()[:maxcol]) for line in content[headrow+1:]],'d') # sort with x as fastest and z as slowest index
N = len(values)
tempGrid = [{},{},{}]
for j in xrange(3):
for i in xrange(N):
tempGrid[j][str(values[i,locol+j])] = True
grid = np.array([len(tempGrid[0]),\
len(tempGrid[1]),\
len(tempGrid[2]),],'i')
dim = np.ones(3)
for i,r in enumerate(grid):
if r > 1:
dim[i] = (max(map(float,tempGrid[i].keys()))-min(map(float,tempGrid[i].keys())))*r/(r-1.0)
if grid[2]==1: # for 2D case set undefined dimension to given unitlength or alternatively give it the length of the smallest element
if options.unitlength == 0.0:
dim[2] = min(dim/grid)
else:
dim[2] = options.unitlength
print dim
if options.undeformed:
Favg = np.eye(3)
else:
Favg = damask.core.math.tensorAvg(
np.reshape(np.transpose(values[:,column['tensor'][options.defgrad]:
column['tensor'][options.defgrad]+9]),
(3,3,grid[0],grid[1],grid[2])))
F = np.reshape(np.transpose(values[:,column['tensor'][options.defgrad]:
column['tensor'][options.defgrad]+9]),
(3,3,grid[0],grid[1],grid[2]))
centroids = damask.core.mesh.deformedCoordsFFT(dim,F,Favg,options.scaling)
nodes = damask.core.mesh.nodesAroundCentres(dim,Favg,centroids)
fields = {\
'tensor': {},\
'vector': {},\
'scalar': {},\
'double': {},\
'triple': {},\
'quadruple': {},\
}
reshape = {\
'tensor': [3,3],\
'vector': [3],\
'scalar': [],\
'double': [2],\
'triple': [3],\
'quadruple': [4],\
}
length = {\
'tensor': 9,\
'vector': 3,\
'scalar': 1,\
'double': 2,\
'triple': 3,\
'quadruple': 4,\
}
# vtk lib out
if False:
points = vtk.vtkPoints()
for z in range (grid[2]+1):
for y in range (grid[1]+1):
for x in range (grid[0]+1):
points.InsertNextPoint(nodes[:,x,y,z])
data=[]
j=0
for datatype in fields.keys():
for what in eval('options.'+datatype):
for label in matches[datatype][what]:
col = column[datatype][label]
if col != -1:
data.append(vtk.vtkFloatArray())
data[j].SetNumberOfComponents(length[datatype])
for i in xrange(grid[2]*grid[1]*grid[0]):
for k in xrange(length[datatype]):
data[j].InsertNextValue(values2[i,col+k])
data[j].SetName(label)
j+=1
if options.output_mesh:
hexs = vtk.vtkCellArray()
i = 0
elems=[]
for z in range (grid[2]):
for y in range (grid[1]):
for x in range (grid[0]):
elems.append(vtk.vtkHexahedron())
base = z*(grid[1]+1)*(grid[0]+1)+y*(grid[0]+1)+x
elems[i].GetPointIds().SetId(0, base)
elems[i].GetPointIds().SetId(1, base+1)
elems[i].GetPointIds().SetId(2, base+grid[0]+2)
elems[i].GetPointIds().SetId(3, base+grid[0]+1)
elems[i].GetPointIds().SetId(4, base+(grid[1]+1)*(grid[0]+1))
elems[i].GetPointIds().SetId(5, base+(grid[1]+1)*(grid[0]+1)+1)
elems[i].GetPointIds().SetId(6, base+(grid[1]+1)*(grid[0]+1)+grid[0]+2)
elems[i].GetPointIds().SetId(7, base+(grid[1]+1)*(grid[0]+1)+grid[0]+1)
hexs.InsertNextCell(elems[i])
i+=1
uGrid = vtk.vtkUnstructuredGrid()
uGrid.SetPoints(points)
i = 0
for z in range (grid[2]):
for y in range (grid[1]):
for x in range (grid[0]):
uGrid.InsertNextCell(elems[i].GetCellType(), elems[i].GetPointIds())
i+=1
for i in xrange(len(data)):
uGrid.GetCellData().AddArray(data[i])
outWriter = vtk.vtkXMLUnstructuredGridWriter()
outWriter.SetDataModeToBinary()
outWriter.SetCompressorTypeToZLib()
(head,tail) = os.path.split(filename)
outWriter.SetFileName(os.path.join(head,'mesh_'+os.path.splitext(tail)[0]+'.vtu'))
outWriter.SetInput(uGrid)
outWriter.Write()
for datatype in fields.keys():
print '\n%s:'%datatype,
fields[datatype]['_order_'] = []
for what in eval('options.'+datatype):
for label in matches[datatype][what]:
col = column[datatype][label]
if col != -1:
print label,
fields[datatype][label] = np.reshape(values[:,col:col+length[datatype]],[grid[0],grid[1],grid[2]]+reshape[datatype])
fields[datatype]['_order_'] += [label]
print '\n'
out = {}
if options.output_mesh: out['mesh'] = vtk_writeASCII_mesh(nodes,fields,grid,sep[options.separator])
if options.output_points: out['points'] = vtk_writeASCII_points(centroids,fields,grid,sep[options.separator])
for what in out.keys():
print what
(head,tail) = os.path.split(filename)
vtk = open(os.path.join(head,what+'_'+os.path.splitext(tail)[0]+'.vtk'), 'w')
output(out[what],{'filepointer':vtk},'File')
vtk.close()
print

32
processing/post/addCalculation.py
View File

@ -74,32 +74,30 @@ for name in filenames:
}
# ------------------------------------------ Evaluate condition ---------------------------------------
if options.condition:
if options.condition is not None:
interpolator = []
condition = options.condition # copy per file, since might be altered inline
condition = options.condition # copy per file, since might be altered inline
breaker = False
for position,operand in enumerate(set(re.findall(r'#(([s]#)?(.+?))#',condition))): # find three groups
for position,operand in enumerate(set(re.findall(r'#(([s]#)?(.+?))#',condition))): # find three groups
condition = condition.replace('#'+operand[0]+'#',
{ '': '{%i}'%position,
's#':'"{%i}"'%position}[operand[1]])
if operand[2] in specials: # special label
if operand[2] in specials: # special label
interpolator += ['specials["%s"]'%operand[2]]
else:
try:
interpolator += ['%s(table.data[%i])'%({ '':'float',
's#':'str'}[operand[1]],
table.label_index(operand[2]))] # ccould be generalized to indexrange as array lookup
table.label_index(operand[2]))] # could be generalized to indexrange as array lookup
except:
damask.util.croak('column "{}" not found.'.format(operand[2]))
breaker = True
if breaker: continue # found mistake in condition evaluation --> next file
if breaker: continue # found mistake in condition evaluation --> next file
evaluator_condition = "'" + condition + "'.format(" + ','.join(interpolator) + ")"
else: condition = ''
# ------------------------------------------ build formulae ----------------------------------------
evaluator = {}
@ -165,19 +163,19 @@ for name in filenames:
for label in output.labels():
oldIndices = table.label_indexrange(label)
Nold = max(1,len(oldIndices)) # Nold could be zero for new columns
Nold = max(1,len(oldIndices)) # Nold could be zero for new columns
Nnew = len(output.label_indexrange(label))
output.data_append(eval(evaluator[label]) if label in options.labels and
(condition == '' or eval(eval(evaluator_condition)))
else np.tile([table.data[i] for i in oldIndices]
if label in tabLabels
else np.nan,
np.ceil(float(Nnew)/Nold))[:Nnew]) # spread formula result into given number of columns
(options.condition is None or eval(eval(evaluator_condition)))
else np.tile([table.data[i] for i in oldIndices]
if label in tabLabels
else np.nan,
np.ceil(float(Nnew)/Nold))[:Nnew]) # spread formula result into given number of columns
outputAlive = output.data_write() # output processed line
outputAlive = output.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------
table.input_close() # close ASCII tables
output.close() # close ASCII tables
table.input_close() # close ASCII tables
output.close() # close ASCII tables

305
processing/post/addCompatibilityMismatch.py
View File

@ -1,14 +1,210 @@
#!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*-
import os,sys
import os
import math
import numpy as np
import scipy.ndimage
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
def cell2node(cellData,grid):
nodeData = 0.0
datalen = np.array(cellData.shape[3:]).prod()
for i in xrange(datalen):
node = scipy.ndimage.convolve(cellData.reshape(tuple(grid[::-1])+(datalen,))[...,i],
np.ones((2,2,2))/8., # 2x2x2 neighborhood of cells
mode = 'wrap',
origin = -1, # offset to have cell origin as center
) # now averaged at cell origins
node = np.append(node,node[np.newaxis,0,:,:,...],axis=0) # wrap along z
node = np.append(node,node[:,0,np.newaxis,:,...],axis=1) # wrap along y
node = np.append(node,node[:,:,0,np.newaxis,...],axis=2) # wrap along x
nodeData = node[...,np.newaxis] if i==0 else np.concatenate((nodeData,node[...,np.newaxis]),axis=-1)
return nodeData
#--------------------------------------------------------------------------------------------------
def deformationAvgFFT(F,grid,size,nodal=False,transformed=False):
"""calculate average cell center (or nodal) deformation for deformation gradient field specified in each grid cell"""
if nodal:
x, y, z = np.meshgrid(np.linspace(0,size[2],1+grid[2]),
np.linspace(0,size[1],1+grid[1]),
np.linspace(0,size[0],1+grid[0]),
indexing = 'ij')
else:
x, y, z = np.meshgrid(np.linspace(size[2]/grid[2]/2.,size[2]-size[2]/grid[2]/2.,grid[2]),
np.linspace(size[1]/grid[1]/2.,size[1]-size[1]/grid[1]/2.,grid[1]),
np.linspace(size[0]/grid[0]/2.,size[0]-size[0]/grid[0]/2.,grid[0]),
indexing = 'ij')
origCoords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
F_fourier = F if transformed else np.fft.rfftn(F,axes=(0,1,2)) # transform or use provided data
Favg = np.real(F_fourier[0,0,0,:,:])/grid.prod() # take zero freq for average
avgDeformation = np.einsum('ml,ijkl->ijkm',Favg,origCoords) # dX = Favg.X
return avgDeformation
#--------------------------------------------------------------------------------------------------
def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
"""calculate cell center (or nodal) displacement for deformation gradient field specified in each grid cell"""
integrator = 0.5j * size / math.pi
kk, kj, ki = np.meshgrid(np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2])),
np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1])),
np.arange(grid[0]//2+1),
indexing = 'ij')
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3)
k_sSquared = np.einsum('...l,...l',k_s,k_s)
k_sSquared[0,0,0] = 1.0 # ignore global average frequency
#--------------------------------------------------------------------------------------------------
# integration in Fourier space
displacement_fourier = -np.einsum('ijkml,ijkl,l->ijkm',
F if transformed else np.fft.rfftn(F,axes=(0,1,2)),
k_s,
integrator,
) / k_sSquared[...,np.newaxis]
#--------------------------------------------------------------------------------------------------
# backtransformation to real space
displacement = np.fft.irfftn(displacement_fourier,grid[::-1],axes=(0,1,2))
return cell2node(displacement,grid) if nodal else displacement
def volTetrahedron(coords):
"""
Return the volume of the tetrahedron with given vertices or sides.
Ifvertices are given they must be in a NumPy array with shape (4,3): the
position vectors of the 4 vertices in 3 dimensions; if the six sides are
given, they must be an array of length 6. If both are given, the sides
will be used in the calculation.
This method implements
Tartaglia's formula using the Cayley-Menger determinant:
|0 1 1 1 1 |
|1 0 s1^2 s2^2 s3^2|
288 V^2 = |1 s1^2 0 s4^2 s5^2|
|1 s2^2 s4^2 0 s6^2|
|1 s3^2 s5^2 s6^2 0 |
where s1, s2, ..., s6 are the tetrahedron side lengths.
from http://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron
"""
# The indexes of rows in the vertices array corresponding to all
# possible pairs of vertices
vertex_pair_indexes = np.array(((0, 1), (0, 2), (0, 3),
(1, 2), (1, 3), (2, 3)))
# Get all the squares of all side lengths from the differences between
# the 6 different pairs of vertex positions
vertices = np.concatenate((coords[0],coords[1],coords[2],coords[3])).reshape([4,3])
vertex1, vertex2 = vertex_pair_indexes[:,0], vertex_pair_indexes[:,1]
sides_squared = np.sum((vertices[vertex1] - vertices[vertex2])**2,axis=-1)
# Set up the Cayley-Menger determinant
M = np.zeros((5,5))
# Fill in the upper triangle of the matrix
M[0,1:] = 1
# The squared-side length elements can be indexed using the vertex
# pair indices (compare with the determinant illustrated above)
M[tuple(zip(*(vertex_pair_indexes + 1)))] = sides_squared
# The matrix is symmetric, so we can fill in the lower triangle by
# adding the transpose
M = M + M.T
return np.sqrt(np.linalg.det(M) / 288)
def volumeMismatch(size,F,nodes):
"""
calculates the volume mismatch
volume mismatch is defined as the difference between volume of reconstructed
(compatible) cube and determinant of defgrad at the FP
"""
coords = np.empty([8,3])
vMismatch = np.empty(grid[::-1])
volInitial = size.prod()/grid.prod()
#--------------------------------------------------------------------------------------------------
# calculate actual volume and volume resulting from deformation gradient
for k in xrange(grid[2]):
for j in xrange(grid[1]):
for i in xrange(grid[0]):
coords[0,0:3] = nodes[k, j, i ,0:3]
coords[1,0:3] = nodes[k ,j, i+1,0:3]
coords[2,0:3] = nodes[k ,j+1,i+1,0:3]
coords[3,0:3] = nodes[k, j+1,i ,0:3]
coords[4,0:3] = nodes[k+1,j, i ,0:3]
coords[5,0:3] = nodes[k+1,j, i+1,0:3]
coords[6,0:3] = nodes[k+1,j+1,i+1,0:3]
coords[7,0:3] = nodes[k+1,j+1,i ,0:3]
vMismatch[k,j,i] = \
( abs(volTetrahedron([coords[6,0:3],coords[0,0:3],coords[7,0:3],coords[3,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[0,0:3],coords[7,0:3],coords[4,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[0,0:3],coords[2,0:3],coords[3,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[0,0:3],coords[2,0:3],coords[1,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[5,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[0,0:3]]))) \
/np.linalg.det(F[k,j,i,0:3,0:3])
return vMismatch/volInitial
def shapeMismatch(size,F,nodes,centres):
"""
Routine to calculate the shape mismatch
shape mismatch is defined as difference between the vectors from the central point to
the corners of reconstructed (combatible) volume element and the vectors calculated by deforming
the initial volume element with the current deformation gradient
"""
coordsInitial = np.empty([8,3])
sMismatch = np.empty(grid[::-1])
#--------------------------------------------------------------------------------------------------
# initial positions
coordsInitial[0,0:3] = [-size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[1,0:3] = [+size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[2,0:3] = [+size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[3,0:3] = [-size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[4,0:3] = [-size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]]
coordsInitial[5,0:3] = [+size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]]
coordsInitial[6,0:3] = [+size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
coordsInitial[7,0:3] = [-size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
coordsInitial = coordsInitial/2.0
#--------------------------------------------------------------------------------------------------
# compare deformed original and deformed positions to actual positions
for k in xrange(grid[2]):
for j in xrange(grid[1]):
for i in xrange(grid[0]):
sMismatch[k,j,i] = \
+ np.linalg.norm(nodes[k, j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[0,0:3]))\
+ np.linalg.norm(nodes[k, j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[1,0:3]))\
+ np.linalg.norm(nodes[k, j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[2,0:3]))\
+ np.linalg.norm(nodes[k, j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[3,0:3]))\
+ np.linalg.norm(nodes[k+1,j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[4,0:3]))\
+ np.linalg.norm(nodes[k+1,j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[5,0:3]))\
+ np.linalg.norm(nodes[k+1,j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[6,0:3]))\
+ np.linalg.norm(nodes[k+1,j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[7,0:3]))
return sMismatch
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
@ -64,79 +260,78 @@ for name in filenames:
errors = []
remarks = []
if table.label_dimension(options.pos) != 3:
errors.append('coordinates "{}" are not a vector.'.format(options.pos))
else: colCoord = table.label_index(options.pos)
if table.label_dimension(options.defgrad) != 9:
errors.append('deformation gradient "{}" is not a tensor.'.format(options.defgrad))
else: colF = table.label_index(options.defgrad)
errors.append('deformation gradient "{}" is not a 3x3 tensor.'.format(options.defgrad))
coordDim = table.label_dimension(options.pos)
if not 3 >= coordDim >= 1:
errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
elif coordDim < 3:
remarks.append('appending {} dimension{} to coordinates "{}"...'.format(3-coordDim,
's' if coordDim < 2 else '',
options.pos))
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
table.close(dismiss=True)
continue
# ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
if options.shape: table.labels_append('shapeMismatch({})'.format(options.defgrad))
if options.volume: table.labels_append('volMismatch({})'.format(options.defgrad))
# --------------- figure out size and grid ---------------------------------------------------------
table.data_readArray()
table.data_readArray([options.defgrad,options.pos])
table.data_rewind()
coords = [np.unique(table.data[:,colCoord+i]) for i in xrange(3)]
if len(table.data.shape) < 2: table.data.shape += (1,) # expand to 2D shape
if table.data[:,9:].shape[1] < 3:
table.data = np.hstack((table.data,
np.zeros((table.data.shape[0],
3-table.data[:,9:].shape[1]),dtype='f'))) # fill coords up to 3D with zeros
coords = [np.unique(table.data[:,9+i]) for i in xrange(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # grid==1 spacing set to smallest among other ones
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings
N = grid.prod()
# --------------- figure out columns to process ---------------------------------------------------
key = '1_'+options.defgrad
if table.label_index(key) == -1:
damask.util.croak('column "{}" not found...'.format(key))
if N != len(table.data): errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*grid))
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
else:
column = table.label_index(key) # remember columns of requested data
# ------------------------------------------ assemble header ---------------------------------------
if options.shape: table.labels_append(['shapeMismatch({})'.format(options.defgrad)])
if options.volume: table.labels_append(['volMismatch({})'.format(options.defgrad)])
table.head_write()
# ------------------------------------------ read deformation gradient field -----------------------
table.data_rewind()
F = np.zeros(N*9,'d').reshape([3,3]+list(grid))
idx = 0
while table.data_read():
(x,y,z) = damask.util.gridLocation(idx,grid) # figure out (x,y,z) position from line count
idx += 1
F[0:3,0:3,x,y,z] = np.array(map(float,table.data[column:column+9]),'d').reshape(3,3)
Favg = damask.core.math.tensorAvg(F)
centres = damask.core.mesh.deformedCoordsFFT(size,F,Favg,[1.0,1.0,1.0])
nodes = damask.core.mesh.nodesAroundCentres(size,Favg,centres)
if options.shape: shapeMismatch = damask.core.mesh.shapeMismatch( size,F,nodes,centres)
if options.volume: volumeMismatch = damask.core.mesh.volumeMismatch(size,F,nodes)
# -----------------------------process data and assemble header -------------------------------------
# ------------------------------------------ process data ------------------------------------------
table.data_rewind()
idx = 0
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
(x,y,z) = damask.util.gridLocation(idx,grid) # figure out (x,y,z) position from line count
idx += 1
if options.shape: table.data_append( shapeMismatch[x,y,z])
if options.volume: table.data_append(volumeMismatch[x,y,z])
outputAlive = table.data_write()
F_fourier = np.fft.rfftn(table.data[:,:9].reshape(grid[2],grid[1],grid[0],3,3),axes=(0,1,2)) # perform transform only once...
nodes = displacementFluctFFT(F_fourier,grid,size,True,transformed=True)\
+ deformationAvgFFT (F_fourier,grid,size,True,transformed=True)
# ------------------------------------------ output finalization -----------------------------------
if options.shape:
table.labels_append(['shapeMismatch({})'.format(options.defgrad)])
centres = displacementFluctFFT(F_fourier,grid,size,False,transformed=True)\
+ deformationAvgFFT (F_fourier,grid,size,False,transformed=True)
if options.volume:
table.labels_append(['volMismatch({})'.format(options.defgrad)])
table.head_write()
if options.shape:
shapeMismatch = shapeMismatch( size,table.data[:,:9].reshape(grid[2],grid[1],grid[0],3,3),nodes,centres)
if options.volume:
volumeMismatch = volumeMismatch(size,table.data[:,:9].reshape(grid[2],grid[1],grid[0],3,3),nodes)
# ------------------------------------------ output data -------------------------------------------
for i in xrange(grid[2]):
for j in xrange(grid[1]):
for k in xrange(grid[0]):
table.data_read()
if options.shape: table.data_append(shapeMismatch[i,j,k])
if options.volume: table.data_append(volumeMismatch[i,j,k])
table.data_write()
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables

6
processing/post/addCurl.py
View File

@ -18,7 +18,7 @@ def curlFFT(geomdim,field):
if n == 3: dataType = 'vector'
elif n == 9: dataType = 'tensor'
field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2),s=shapeFFT)
field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
curl_fourier = np.empty(field_fourier.shape,'c16')
# differentiation in Fourier space
@ -56,7 +56,7 @@ def curlFFT(geomdim,field):
curl_fourier[i,j,k,2] = ( field_fourier[i,j,k,1]*xi[0]\
-field_fourier[i,j,k,0]*xi[1]) *TWOPIIMG
return np.fft.fftpack.irfftn(curl_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n])
return np.fft.irfftn(curl_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n])
# --------------------------------------------------------------------
@ -158,7 +158,7 @@ for name in filenames:
# we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation
stack.append(curlFFT(size[::-1],
table.data[:,data['column'][i]:data['column'][i]+data['dim']].
reshape([grid[2],grid[1],grid[0]]+data['shape'])))
reshape(grid[::-1].tolist()+data['shape'])))
# ------------------------------------------ output result -----------------------------------------

34
processing/post/addDisplacement.py
View File

@ -18,7 +18,7 @@ def cell2node(cellData,grid):
datalen = np.array(cellData.shape[3:]).prod()
for i in xrange(datalen):
node = scipy.ndimage.convolve(cellData.reshape(tuple(grid)+(datalen,))[...,i],
node = scipy.ndimage.convolve(cellData.reshape(tuple(grid[::-1])+(datalen,))[...,i],
np.ones((2,2,2))/8., # 2x2x2 neighborhood of cells
mode = 'wrap',
origin = -1, # offset to have cell origin as center
@ -35,14 +35,14 @@ def cell2node(cellData,grid):
def displacementAvgFFT(F,grid,size,nodal=False,transformed=False):
"""calculate average cell center (or nodal) displacement for deformation gradient field specified in each grid cell"""
if nodal:
x, y, z = np.meshgrid(np.linspace(0,size[0],1+grid[0]),
x, y, z = np.meshgrid(np.linspace(0,size[2],1+grid[2]),
np.linspace(0,size[1],1+grid[1]),
np.linspace(0,size[2],1+grid[2]),
np.linspace(0,size[0],1+grid[0]),
indexing = 'ij')
else:
x, y, z = np.meshgrid(np.linspace(0,size[0],grid[0],endpoint=False),
x, y, z = np.meshgrid(np.linspace(0,size[2],grid[2],endpoint=False),
np.linspace(0,size[1],grid[1],endpoint=False),
np.linspace(0,size[2],grid[2],endpoint=False),
np.linspace(0,size[0],grid[0],endpoint=False),
indexing = 'ij')
origCoords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
@ -69,7 +69,7 @@ def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
#--------------------------------------------------------------------------------------------------
# integration in Fourier space
displacement_fourier = +np.einsum('ijkml,ijkl,l->ijkm',
displacement_fourier = -np.einsum('ijkml,ijkl,l->ijkm',
F if transformed else np.fft.rfftn(F,axes=(0,1,2)),
k_s,
integrator,
@ -78,7 +78,7 @@ def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
#--------------------------------------------------------------------------------------------------
# backtransformation to real space
displacement = np.fft.irfftn(displacement_fourier,grid,axes=(0,1,2))
displacement = np.fft.irfftn(displacement_fourier,grid[::-1],axes=(0,1,2))
return cell2node(displacement,grid) if nodal else displacement
@ -186,8 +186,8 @@ for name in filenames:
F_fourier = np.fft.rfftn(table.data[:,:9].reshape(grid[2],grid[1],grid[0],3,3),axes=(0,1,2)) # perform transform only once...
displacement = displacementFluctFFT(F_fourier,grid,size,options.nodal,transformed=True)
avgDisplacement = displacementAvgFFT (F_fourier,grid,size,options.nodal,transformed=True)
fluctDisplacement = displacementFluctFFT(F_fourier,grid,size,options.nodal,transformed=True)
avgDisplacement = displacementAvgFFT (F_fourier,grid,size,options.nodal,transformed=True)
# ------------------------------------------ assemble header ---------------------------------------
@ -203,18 +203,18 @@ for name in filenames:
# ------------------------------------------ output data -------------------------------------------
zrange = np.linspace(0,size[2],1+grid[2]) if options.nodal else xrange(grid[2])
yrange = np.linspace(0,size[1],1+grid[1]) if options.nodal else xrange(grid[1])
xrange = np.linspace(0,size[0],1+grid[0]) if options.nodal else xrange(grid[0])
Zrange = np.linspace(0,size[2],1+grid[2]) if options.nodal else xrange(grid[2])
Yrange = np.linspace(0,size[1],1+grid[1]) if options.nodal else xrange(grid[1])
Xrange = np.linspace(0,size[0],1+grid[0]) if options.nodal else xrange(grid[0])
for i,z in enumerate(zrange):
for j,y in enumerate(yrange):
for k,x in enumerate(xrange):
for i,z in enumerate(Zrange):
for j,y in enumerate(Yrange):
for k,x in enumerate(Xrange):
if options.nodal: table.data_clear()
else: table.data_read()
table.data_append([x,y,z] if options.nodal else [])
table.data_append(list(avgDisplacement[i,j,k,:]))
table.data_append(list( displacement[i,j,k,:]))
table.data_append(list( avgDisplacement[i,j,k,:]))
table.data_append(list(fluctDisplacement[i,j,k,:]))
table.data_write()
# ------------------------------------------ output finalization -----------------------------------

6
processing/post/addDivergence.py
View File

@ -15,7 +15,7 @@ def divFFT(geomdim,field):
N = grid.prod() # field size
n = np.array(np.shape(field)[3:]).prod() # data size
field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2),s=shapeFFT)
field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
div_fourier = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16') # size depents on whether tensor or vector
# differentiation in Fourier space
@ -42,7 +42,7 @@ def divFFT(geomdim,field):
elif n == 3: # vector, 3 -> 1
div_fourier[i,j,k] = sum(field_fourier[i,j,k,0:3]*xi) *TWOPIIMG
return np.fft.fftpack.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3])
return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3])
# --------------------------------------------------------------------
@ -145,7 +145,7 @@ for name in filenames:
# we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation
stack.append(divFFT(size[::-1],
table.data[:,data['column'][i]:data['column'][i]+data['dim']].
reshape([grid[2],grid[1],grid[0]]+data['shape'])))
reshape(grid[::-1].tolist()+data['shape'])))
# ------------------------------------------ output result -----------------------------------------

6
processing/post/addGradient.py
View File

@ -18,7 +18,7 @@ def gradFFT(geomdim,field):
if n == 3: dataType = 'vector'
elif n == 1: dataType = 'scalar'
field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2),s=shapeFFT)
field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
grad_fourier = np.empty(field_fourier.shape+(3,),'c16')
# differentiation in Fourier space
@ -46,7 +46,7 @@ def gradFFT(geomdim,field):
grad_fourier[i,j,k,1,:] = field_fourier[i,j,k,1]*xi *TWOPIIMG # tensor field from vector data
grad_fourier[i,j,k,2,:] = field_fourier[i,j,k,2]*xi *TWOPIIMG
return np.fft.fftpack.irfftn(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n])
return np.fft.irfftn(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n])
# --------------------------------------------------------------------
@ -148,7 +148,7 @@ for name in filenames:
# we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation
stack.append(gradFFT(size[::-1],
table.data[:,data['column'][i]:data['column'][i]+data['dim']].
reshape([grid[2],grid[1],grid[0]]+data['shape'])))
reshape(grid[::-1].tolist()+data['shape'])))
# ------------------------------------------ output result -----------------------------------------