This commit is contained in:
Martin Diehl 2019-03-09 07:47:01 +01:00
parent 0d08659b2a
commit 66e6a6ec68
7 changed files with 99 additions and 139 deletions

View File

@ -90,7 +90,6 @@ subroutine CPFEM_init
use prec, only: &
pInt, pReal
use IO, only: &
IO_timeStamp, &
IO_error
use numerics, only: &
worldrank

View File

@ -7,11 +7,6 @@
!> results
!--------------------------------------------------------------------------------------------------
program DAMASK_FEM
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
#include <petsc/finclude/petscsys.h>
use PetscDM
use prec, only: &
@ -31,8 +26,7 @@ program DAMASK_FEM
IO_error, &
IO_lc, &
IO_intOut, &
IO_warning, &
IO_timeStamp
IO_warning
use math ! need to include the whole module for FFTW
use CPFEM2, only: &
CPFEM_initAll
@ -118,8 +112,6 @@ program DAMASK_FEM
! init DAMASK (all modules)
call CPFEM_initAll
write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
! reading basic information from load case file and allocate data structure containing load cases
call DMGetDimension(geomMesh,dimPlex,ierr)! CHKERRQ(ierr) !< dimension of mesh (2D or 3D)

View File

@ -66,9 +66,7 @@ contains
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine FEM_mech_init(fieldBC)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use IO, only: &
IO_timeStamp, &
IO_error
use DAMASK_interface, only: &
getSolverJobName
@ -111,8 +109,6 @@ subroutine FEM_mech_init(fieldBC)
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
! Setup FEM mech mesh

View File

@ -36,17 +36,11 @@ contains
!> @brief allocates all neccessary fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
subroutine homogenization_isostrain_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use debug, only: &
debug_HOMOGENIZATION, &
debug_level, &
debug_levelBasic
use IO, only: &
IO_timeStamp, &
IO_error
use material, only: &
homogenization_type, &
@ -67,8 +61,6 @@ subroutine homogenization_isostrain_init()
tag = ''
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_ISOSTRAIN_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
Ninstance = int(count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID),pInt)
if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) &

View File

@ -6,11 +6,11 @@
!--------------------------------------------------------------------------------------------------
module homogenization_none
implicit none
private
public :: &
homogenization_none_init
implicit none
private
public :: &
homogenization_none_init
contains
@ -18,52 +18,42 @@ contains
!> @brief allocates all neccessary fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
subroutine homogenization_none_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use prec, only: &
pInt
use debug, only: &
debug_HOMOGENIZATION, &
debug_level, &
debug_levelBasic
use IO, only: &
IO_timeStamp
use material, only: &
homogenization_type, &
material_homog, &
homogState, &
HOMOGENIZATION_NONE_LABEL, &
HOMOGENIZATION_NONE_ID
implicit none
integer(pInt) :: &
Ninstance, &
h, &
NofMyHomog
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_NONE_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
Ninstance = int(count(homogenization_type == HOMOGENIZATION_NONE_ID),pInt)
if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
do h = 1_pInt, size(homogenization_type)
if (homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle
NofMyHomog = count(material_homog == h)
homogState(h)%sizeState = 0_pInt
homogState(h)%sizePostResults = 0_pInt
allocate(homogState(h)%state0 (0_pInt,NofMyHomog))
allocate(homogState(h)%subState0(0_pInt,NofMyHomog))
allocate(homogState(h)%state (0_pInt,NofMyHomog))
enddo
use debug, only: &
debug_HOMOGENIZATION, &
debug_level, &
debug_levelBasic
use config, only: &
config_homogenization
use material, only: &
homogenization_type, &
material_homog, &
homogState, &
HOMOGENIZATION_NONE_LABEL, &
HOMOGENIZATION_NONE_ID
implicit none
integer :: &
Ninstance, &
h, &
NofMyHomog
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_NONE_label//' init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_NONE_ID)
if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
do h = 1, size(homogenization_type)
if (homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle
NofMyHomog = count(material_homog == h)
homogState(h)%sizeState = 0
homogState(h)%sizePostResults = 0
allocate(homogState(h)%state0 (0,NofMyHomog))
allocate(homogState(h)%subState0(0,NofMyHomog))
allocate(homogState(h)%state (0,NofMyHomog))
enddo
end subroutine homogenization_none_init

View File

@ -280,14 +280,8 @@ contains
!> material.config
!--------------------------------------------------------------------------------------------------
subroutine material_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO, only: &
IO_error, &
IO_timeStamp
IO_error
use debug, only: &
debug_level, &
debug_material, &
@ -321,8 +315,6 @@ subroutine material_init()
myDebug = debug_level(debug_material)
write(6,'(/,a)') ' <<<+- material init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
call material_parsePhase()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6)

View File

@ -503,64 +503,63 @@ end subroutine utilities_FFTvectorBackward
!> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierGammaConvolution(fieldAim)
use numerics, only: &
memory_efficient
use math, only: &
math_det33, &
math_invert
use mesh, only: &
grid3, &
grid, &
grid3Offset
implicit none
real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution
complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx
real(pReal) :: matA(6,6), matInvA(6,6)
integer(pInt) :: &
i, j, k, &
l, m, n, o
logical :: err
write(6,'(/,a)') ' ... doing gamma convolution ...............................................'
flush(6)
use numerics, only: &
memory_efficient
use math, only: &
math_det33, &
math_invert2
use mesh, only: &
grid3, &
grid, &
grid3Offset
implicit none
real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution
complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx
real(pReal), dimension(6,6) :: A, A_inv
integer :: &
i, j, k, &
l, m, n, o
logical :: err
write(6,'(/,a)') ' ... doing gamma convolution ...............................................'
flush(6)
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium)
memoryEfficient: if(memory_efficient) then
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red
if (any([i,j,k+grid3Offset] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex)
matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex)
if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then
call math_invert(6_pInt, matA, matInvA, err)
temp33_complex = cmplx(matInvA(1:3,1:3),matInvA(1:3,4:6),pReal)
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt) &
gamma_hat(l,m,n,o,1,1,1) = temp33_complex(l,n)*conjg(-xi1st(o,i,j,k))*xi1st(m,i,j,k)
else
gamma_hat(1:3,1:3,1:3,1:3,1,1,1) = cmplx(0.0_pReal,0.0_pReal,pReal)
endif
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
endif
enddo; enddo; enddo
else memoryEfficient
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorField_fourier(1:3,1:3,i,j,k))
tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
enddo; enddo; enddo
endif memoryEfficient
if (grid3Offset == 0_pInt) &
tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal)
memoryEfficient: if(memory_efficient) then
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
if (any([i,j,k+grid3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall(l = 1:3, m = 1:3) &
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
forall(l = 1:3, m = 1:3) &
temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
A(1:3,1:3) = real(temp33_complex); A(4:6,4:6) = real(temp33_complex)
A(1:3,4:6) = aimag(temp33_complex); A(4:6,1:3) = -aimag(temp33_complex)
if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
call math_invert2(A_inv, err, A)
temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
forall(l=1:3, m=1:3, n=1:3, o=1:3) &
gamma_hat(l,m,n,o,1,1,1) = temp33_complex(l,n)*conjg(-xi1st(o,i,j,k))*xi1st(m,i,j,k)
else
gamma_hat(1:3,1:3,1:3,1:3,1,1,1) = cmplx(0.0_pReal,0.0_pReal,pReal)
endif
forall(l = 1:3, m = 1:3) &
temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
endif
enddo; enddo; enddo
else memoryEfficient
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
forall(l = 1:3, m = 1:3) &
temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorField_fourier(1:3,1:3,i,j,k))
tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
enddo; enddo; enddo
endif memoryEfficient
if (grid3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal)
end subroutine utilities_fourierGammaConvolution
@ -725,7 +724,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
math_99to3333, &
math_rotate_forward3333, &
math_rotate_forward33, &
math_invert
math_invert2
implicit none
real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance
@ -768,7 +767,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
c_reduced(k,j) = temp99_Real(n,m)
endif; enddo; endif; enddo
call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness
call math_invert2(s_reduced, errmatinv, c_reduced) ! invert reduced stiffness
if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
if (errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
temp99_Real = 0.0_pReal ! fill up compliance with zeros