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

Marc precision handling

This commit is contained in:
Martin Diehl 2022-06-27 08:37:41 +00:00 committed by Daniel Otto de Mentock
parent 0e4a8f284d
commit 3b57934d6e
19 changed files with 108 additions and 131 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@ -108,7 +108,7 @@ file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PE
string(REPLACE "PETSC_FC_INCLUDES = " "" PETSC_INCLUDES "${PETSC_INCLUDES}") string(REPLACE "PETSC_FC_INCLUDES = " "" PETSC_INCLUDES "${PETSC_INCLUDES}")
message("PETSC_INCLUDES:\n${PETSC_INCLUDES}\n") message("PETSC_INCLUDES:\n${PETSC_INCLUDES}\n")
set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}") set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS}")
set(CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${CMAKE_Fortran_COMPILER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}") set(CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${CMAKE_Fortran_COMPILER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}")
if(CMAKE_BUILD_TYPE STREQUAL "DEBUG") if(CMAKE_BUILD_TYPE STREQUAL "DEBUG")

7
cmake/Compiler-GNU.cmake
View File

@ -135,10 +135,3 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fsanitize=undefined")
# detect undefined behavior # detect undefined behavior
# Additional options # Additional options
# -fsanitize=address,leak,thread # -fsanitize=address,leak,thread
#------------------------------------------------------------------------------------------------
# precision settings
set (PRECISION_FLAGS "${PRECISION_FLAGS} -fdefault-real-8")
# set precision to 8 bytes for standard real (=8 for pReal). Will set size of double to 16 bytes as long as -fdefault-double-8 is not set
set (PRECISION_FLAGS "${PRECISION_FLAGS} -fdefault-double-8")
# set precision to 8 bytes for double real, would be 16 bytes if -fdefault-real-8 is used

5
cmake/Compiler-Intel.cmake
View File

@ -118,8 +118,3 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all")
# -check: Checks at runtime, where # -check: Checks at runtime, where
# arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays (performance?) # arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays (performance?)
# stack: # stack:
#------------------------------------------------------------------------------------------------
# precision settings
set (PRECISION_FLAGS "${PRECISION_FLAGS} -real-size 64")
# set precision for standard real to 32 | 64 | 128 (= 4 | 8 | 16 bytes, type pReal is always 8 bytes)

5
cmake/Compiler-IntelLLVM.cmake
View File

@ -117,8 +117,3 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all")
# -check: Checks at runtime, where # -check: Checks at runtime, where
# arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays (performance?) # arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays (performance?)
# stack: # stack:
#------------------------------------------------------------------------------------------------
# precision settings
set (PRECISION_FLAGS "${PRECISION_FLAGS} -real-size 64")
# set precision for standard real to 32 | 64 | 128 (= 4 | 8 | 16 bytes, type pReal is always 8 bytes)

24
src/Marc/DAMASK_Marc.f90
View File

@ -211,7 +211,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
use OMP_LIB use OMP_LIB
implicit none implicit none
integer, intent(in) :: & ! according to MSC.Marc 2012 Manual D integer(pI64), intent(in) :: & ! according to MSC.Marc 2012 Manual D
ngens, & !< size of stress-strain law ngens, & !< size of stress-strain law
nn, & !< integration point number nn, & !< integration point number
ndi, & !< number of direct components ndi, & !< number of direct components
@ -224,7 +224,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
jtype, & !< element type jtype, & !< element type
ifr, & !< set to 1 if R has been calculated ifr, & !< set to 1 if R has been calculated
ifu !< set to 1 if stretch has been calculated ifu !< set to 1 if stretch has been calculated
integer, dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D integer(pI64), dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D
m, & !< (1) user element number, (2) internal element number m, & !< (1) user element number, (2) internal element number
matus, & !< (1) user material identification number, (2) internal material identification number matus, & !< (1) user material identification number, (2) internal material identification number
kcus, & !< (1) layer number, (2) internal layer number kcus, & !< (1) layer number, (2) internal layer number
@ -362,7 +362,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
endif endif
lastLovl = lovl lastLovl = lovl
call materialpoint_general(computationMode,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) call materialpoint_general(computationMode,ffn,ffn1,t(1),timinc,int(m(1)),int(nn),stress,ddsdde)
d = ddsdde(1:ngens,1:ngens) d = ddsdde(1:ngens,1:ngens)
s = stress(1:ndi+nshear) s = stress(1:ndi+nshear)
@ -383,16 +383,17 @@ subroutine flux(f,ts,n,time)
use discretization_Marc use discretization_Marc
implicit none implicit none
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
ts ts
integer, dimension(10), intent(in) :: & integer(pI64), dimension(10), intent(in) :: &
n n
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
time time
real(pReal), dimension(2), intent(out) :: & real(pReal), dimension(2), intent(out) :: &
f f
f(1) = homogenization_f_T(discretization_Marc_FEM2DAMASK_cell(n(3),n(1)))
f(1) = homogenization_f_T(discretization_Marc_FEM2DAMASK_cell(int(n(3)),int(n(1))))
f(2) = 0.0_pReal f(2) = 0.0_pReal
end subroutine flux end subroutine flux
@ -410,7 +411,8 @@ subroutine uedinc(inc,incsub)
use discretization_Marc use discretization_Marc
implicit none implicit none
integer, intent(in) :: inc, incsub integer(pI64), intent(in) :: inc, incsub
integer :: n, nqncomp, nqdatatype integer :: n, nqncomp, nqdatatype
integer, save :: inc_written integer, save :: inc_written
real(pReal), allocatable, dimension(:,:) :: d_n real(pReal), allocatable, dimension(:,:) :: d_n
@ -427,9 +429,9 @@ subroutine uedinc(inc,incsub)
enddo enddo
call discretization_Marc_UpdateNodeAndIpCoords(d_n) call discretization_Marc_UpdateNodeAndIpCoords(d_n)
call materialpoint_results(inc,cptim) call materialpoint_results(int(inc),cptim)
inc_written = inc inc_written = int(inc)
endif endif
end subroutine uedinc end subroutine uedinc

4
src/grid/DAMASK_grid.f90
View File

@ -272,7 +272,7 @@ program DAMASK_grid
write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'R:',& write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'R:',&
transpose(loadCases(l)%rot%asMatrix()) transpose(loadCases(l)%rot%asMatrix())
if (loadCases(l)%r <= 0.0) errorID = 833 if (loadCases(l)%r <= 0.0_pReal) errorID = 833
if (loadCases(l)%t < 0.0_pReal) errorID = 834 if (loadCases(l)%t < 0.0_pReal) errorID = 834
if (loadCases(l)%N < 1) errorID = 835 if (loadCases(l)%N < 1) errorID = 835
if (loadCases(l)%f_out < 1) errorID = 836 if (loadCases(l)%f_out < 1) errorID = 836
@ -501,7 +501,7 @@ subroutine getMaskedTensor(values,mask,tensor)
integer :: i,j integer :: i,j
values = 0.0 values = 0.0_pReal
do i = 1,3 do i = 1,3
row => tensor%get(i) row => tensor%get(i)
do j = 1,3 do j = 1,3

2
src/grid/VTI.f90
View File

@ -151,7 +151,7 @@ subroutine VTI_readCellsSizeOrigin(cells,geomSize,origin, &
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
fileContent fileContent
character(len=:), allocatable :: dataType, headerType character(len=:), allocatable :: headerType
logical :: inFile, inImage, compressed logical :: inFile, inImage, compressed
integer(pI64) :: & integer(pI64) :: &
startPos, endPos startPos, endPos

6
src/grid/discretization_grid.f90
View File

@ -229,9 +229,9 @@ pure function cellSurfaceArea(geomSize,cells)
real(pReal), dimension(6,1,product(cells)) :: cellSurfaceArea real(pReal), dimension(6,1,product(cells)) :: cellSurfaceArea
cellSurfaceArea(1:2,1,:) = geomSize(2)/real(cells(2)) * geomSize(3)/real(cells(3)) cellSurfaceArea(1:2,1,:) = geomSize(2)/real(cells(2),pReal) * geomSize(3)/real(cells(3),pReal)
cellSurfaceArea(3:4,1,:) = geomSize(3)/real(cells(3)) * geomSize(1)/real(cells(1)) cellSurfaceArea(3:4,1,:) = geomSize(3)/real(cells(3),pReal) * geomSize(1)/real(cells(1),pReal)
cellSurfaceArea(5:6,1,:) = geomSize(1)/real(cells(1)) * geomSize(2)/real(cells(2)) cellSurfaceArea(5:6,1,:) = geomSize(1)/real(cells(1),pReal) * geomSize(2)/real(cells(2),pReal)
end function cellSurfaceArea end function cellSurfaceArea

20
src/grid/grid_mech_FEM.f90
View File

@ -221,14 +221,14 @@ subroutine grid_mechanical_FEM_init
delta = geomSize/real(cells,pReal) ! grid spacing delta = geomSize/real(cells,pReal) ! grid spacing
detJ = product(delta) ! cell volume detJ = product(delta) ! cell volume
BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & BMat = reshape(real([-delta(1)**(-1),-delta(2)**(-1),-delta(3)**(-1), &
1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & delta(1)**(-1),-delta(2)**(-1),-delta(3)**(-1), &
-1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & -delta(1)**(-1), delta(2)**(-1),-delta(3)**(-1), &
1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & delta(1)**(-1), delta(2)**(-1),-delta(3)**(-1), &
-1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & -delta(1)**(-1),-delta(2)**(-1), delta(3)**(-1), &
1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & delta(1)**(-1),-delta(2)**(-1), delta(3)**(-1), &
-1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3), & -delta(1)**(-1), delta(2)**(-1), delta(3)**(-1), &
1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3)],pReal), [3,8])/4.0_pReal ! shape function derivative matrix delta(1)**(-1), delta(2)**(-1), delta(3)**(-1)],pReal), [3,8])/4.0_pReal ! shape function derivative matrix
HGMat = matmul(transpose(HGcomp),HGcomp) & HGMat = matmul(transpose(HGcomp),HGcomp) &
* HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal ! hourglass stabilization matrix * HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal ! hourglass stabilization matrix
@ -656,7 +656,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
MatNullSpace :: matnull MatNullSpace :: matnull
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
BMatFull = 0.0 BMatFull = 0.0_pReal
BMatFull(1:3,1 :8 ) = BMat BMatFull(1:3,1 :8 ) = BMat
BMatFull(4:6,9 :16) = BMat BMatFull(4:6,9 :16) = BMat
BMatFull(7:9,17:24) = BMat BMatFull(7:9,17:24) = BMat
@ -686,7 +686,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
enddo; enddo; enddo enddo; enddo; enddo
row = col row = col
ce = ce + 1 ce = ce + 1
K_ele = 0.0 K_ele = 0.0_pReal
K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + & K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + &
homogenization_dPdF(2,2,2,2,ce) + & homogenization_dPdF(2,2,2,2,ce) + &
homogenization_dPdF(3,3,3,3,ce))/3.0_pReal homogenization_dPdF(3,3,3,3,ce))/3.0_pReal

28
src/grid/spectral_utilities.f90
View File

@ -198,7 +198,7 @@ subroutine spectral_utilities_init()
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
cells1Red = cells(1)/2 + 1 cells1Red = cells(1)/2 + 1
wgt = 1.0/real(product(cells),pReal) wgt = real(product(cells),pReal)**(-1)
num%memory_efficient = num_grid%get_asInt('memory_efficient', defaultVal=1) > 0 ! ToDo: should be logical in YAML file num%memory_efficient = num_grid%get_asInt('memory_efficient', defaultVal=1) > 0 ! ToDo: should be logical in YAML file
num%divergence_correction = num_grid%get_asInt('divergence_correction', defaultVal=2) num%divergence_correction = num_grid%get_asInt('divergence_correction', defaultVal=2)
@ -265,7 +265,7 @@ subroutine spectral_utilities_init()
N = fftw_mpi_local_size_many(3,[cellsFFTW(3),cellsFFTW(2),cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T],& N = fftw_mpi_local_size_many(3,[cellsFFTW(3),cellsFFTW(2),cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T],&
tensorSize,FFTW_MPI_DEFAULT_BLOCK,PETSC_COMM_WORLD,z,devNull) tensorSize,FFTW_MPI_DEFAULT_BLOCK,PETSC_COMM_WORLD,z,devNull)
if (z /= cells3) error stop 'domain decomposition mismatch (tensor)' if (int(z) /= cells3) error stop 'domain decomposition mismatch (tensor)'
tensorField = fftw_alloc_complex(N) tensorField = fftw_alloc_complex(N)
call c_f_pointer(tensorField,tensorField_real, & call c_f_pointer(tensorField,tensorField_real, &
[3_C_INTPTR_T,3_C_INTPTR_T,2_C_INTPTR_T*(cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T),cellsFFTW(2),z]) [3_C_INTPTR_T,3_C_INTPTR_T,2_C_INTPTR_T*(cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T),cellsFFTW(2),z])
@ -274,7 +274,7 @@ subroutine spectral_utilities_init()
N = fftw_mpi_local_size_many(3,[cellsFFTW(3),cellsFFTW(2),cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T],& N = fftw_mpi_local_size_many(3,[cellsFFTW(3),cellsFFTW(2),cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T],&
vectorSize,FFTW_MPI_DEFAULT_BLOCK,PETSC_COMM_WORLD,z,devNull) vectorSize,FFTW_MPI_DEFAULT_BLOCK,PETSC_COMM_WORLD,z,devNull)
if (z /= cells3) error stop 'domain decomposition mismatch (vector)' if (int(z) /= cells3) error stop 'domain decomposition mismatch (vector)'
vectorField = fftw_alloc_complex(N) vectorField = fftw_alloc_complex(N)
call c_f_pointer(vectorField,vectorField_real, & call c_f_pointer(vectorField,vectorField_real, &
[3_C_INTPTR_T,2_C_INTPTR_T*(cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T),cellsFFTW(2),z]) [3_C_INTPTR_T,2_C_INTPTR_T*(cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T),cellsFFTW(2),z])
@ -283,7 +283,7 @@ subroutine spectral_utilities_init()
N = fftw_mpi_local_size_3d(cellsFFTW(3),cellsFFTW(2),cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T,& N = fftw_mpi_local_size_3d(cellsFFTW(3),cellsFFTW(2),cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T,&
PETSC_COMM_WORLD,z,devNull) PETSC_COMM_WORLD,z,devNull)
if (z /= cells3) error stop 'domain decomposition mismatch (scalar)' if (int(z) /= cells3) error stop 'domain decomposition mismatch (scalar)'
scalarField = fftw_alloc_complex(N) scalarField = fftw_alloc_complex(N)
call c_f_pointer(scalarField,scalarField_real, & call c_f_pointer(scalarField,scalarField_real, &
[2_C_INTPTR_T*(cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T),cellsFFTW(2),z]) [2_C_INTPTR_T*(cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T),cellsFFTW(2),z])
@ -385,17 +385,17 @@ subroutine utilities_updateGamma(C)
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
end do end do
do concurrent(l = 1:3, m = 1:3) do concurrent(l = 1:3, m = 1:3)
temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal,pReal)*xiDyad_cmplx)
end do end do
#else #else
forall(l = 1:3, m = 1:3) & forall(l = 1:3, m = 1:3) &
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
forall(l = 1:3, m = 1:3) & forall(l = 1:3, m = 1:3) &
temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal,pReal)*xiDyad_cmplx)
#endif #endif
A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) = temp33_cmplx%re A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) = temp33_cmplx%re
A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im
if (abs(math_det33(A(1:3,1:3))) > 1e-16) then if (abs(math_det33(A(1:3,1:3))) > 1.e-16_pReal) then
call math_invert(A_inv, err, A) call math_invert(A_inv, err, A)
temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal) temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
#ifndef __INTEL_COMPILER #ifndef __INTEL_COMPILER
@ -518,17 +518,17 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
end do end do
do concurrent(l = 1:3, m = 1:3) do concurrent(l = 1:3, m = 1:3)
temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal,pReal)*xiDyad_cmplx)
end do end do
#else #else
forall(l = 1:3, m = 1:3) & forall(l = 1:3, m = 1:3) &
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
forall(l = 1:3, m = 1:3) & forall(l = 1:3, m = 1:3) &
temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal,pReal)*xiDyad_cmplx)
#endif #endif
A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) = temp33_cmplx%re A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) = temp33_cmplx%re
A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im
if (abs(math_det33(A(1:3,1:3))) > 1e-16) then if (abs(math_det33(A(1:3,1:3))) > 1.e-16_pReal) then
call math_invert(A_inv, err, A) call math_invert(A_inv, err, A)
temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal) temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
#ifndef __INTEL_COMPILER #ifndef __INTEL_COMPILER
@ -587,8 +587,8 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t)
!$OMP PARALLEL DO PRIVATE(GreenOp_hat) !$OMP PARALLEL DO PRIVATE(GreenOp_hat)
do k = 1, cells3; do j = 1, cells(2) ;do i = 1, cells1Red do k = 1, cells3; do j = 1, cells(2) ;do i = 1, cells1Red
GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) & GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) &
/ (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal) & / (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal,pReal) &
* sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k)))) * sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal,pReal),xi1st(1:3,i,j,k))))
scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -608,7 +608,7 @@ real(pReal) function utilities_divergenceRMS()
print'(/,1x,a)', '... calculating divergence ................................................' print'(/,1x,a)', '... calculating divergence ................................................'
flush(IO_STDOUT) flush(IO_STDOUT)
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space ! calculating RMS divergence criterion in Fourier space
@ -652,7 +652,7 @@ real(pReal) function utilities_curlRMS()
print'(/,1x,a)', '... calculating curl ......................................................' print'(/,1x,a)', '... calculating curl ......................................................'
flush(IO_STDOUT) flush(IO_STDOUT)
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating max curl criterion in Fourier space ! calculating max curl criterion in Fourier space

6
src/homogenization_mechanical_RGC.f90
View File

@ -601,9 +601,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
! calculate the stress and penalty due to volume discrepancy ! calculate the stress and penalty due to volume discrepancy
vPen = 0.0_pReal vPen = 0.0_pReal
do i = 1,nGrain do i = 1,nGrain
vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* & vPen(:,:,i) = -real(nGrain,pReal)**(-1)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr &
sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* & * sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0_pReal),vDiscrep) &
gVol(i)*transpose(math_inv33(fDef(:,:,i))) * gVol(i)*transpose(math_inv33(fDef(:,:,i)))
end do end do
end subroutine volumePenalty end subroutine volumePenalty

8
src/lattice.f90
View File

@ -2008,7 +2008,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
],pReal),shape(CFTOHP_SYSTEMTRANS)) ],pReal),shape(CFTOHP_SYSTEMTRANS))
real(pReal), dimension(4,cF_Ntrans), parameter :: & real(pReal), dimension(4,cF_Ntrans), parameter :: &
CFTOCI_SYSTEMTRANS = reshape([& CFTOCI_SYSTEMTRANS = real(reshape([&
0.0, 1.0, 0.0, 10.26, & ! Pitsch OR (Ma & Hartmaier 2014, Table 3) 0.0, 1.0, 0.0, 10.26, & ! Pitsch OR (Ma & Hartmaier 2014, Table 3)
0.0,-1.0, 0.0, 10.26, & 0.0,-1.0, 0.0, 10.26, &
0.0, 0.0, 1.0, 10.26, & 0.0, 0.0, 1.0, 10.26, &
@ -2021,7 +2021,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
-1.0, 0.0, 0.0, 10.26, & -1.0, 0.0, 0.0, 10.26, &
0.0, 1.0, 0.0, 10.26, & 0.0, 1.0, 0.0, 10.26, &
0.0,-1.0, 0.0, 10.26 & 0.0,-1.0, 0.0, 10.26 &
],shape(CFTOCI_SYSTEMTRANS)) ],shape(CFTOCI_SYSTEMTRANS)),pReal)
integer, dimension(9,cF_Ntrans), parameter :: & integer, dimension(9,cF_Ntrans), parameter :: &
CFTOCI_BAINVARIANT = reshape( [& CFTOCI_BAINVARIANT = reshape( [&
@ -2040,7 +2040,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
],shape(CFTOCI_BAINVARIANT)) ],shape(CFTOCI_BAINVARIANT))
real(pReal), dimension(4,cF_Ntrans), parameter :: & real(pReal), dimension(4,cF_Ntrans), parameter :: &
CFTOCI_BAINROT = reshape([& CFTOCI_BAINROT = real(reshape([&
1.0, 0.0, 0.0, 45.0, & ! Rotate cF austensite to bain variant 1.0, 0.0, 0.0, 45.0, & ! Rotate cF austensite to bain variant
1.0, 0.0, 0.0, 45.0, & 1.0, 0.0, 0.0, 45.0, &
1.0, 0.0, 0.0, 45.0, & 1.0, 0.0, 0.0, 45.0, &
@ -2053,7 +2053,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
0.0, 0.0, 1.0, 45.0, & 0.0, 0.0, 1.0, 45.0, &
0.0, 0.0, 1.0, 45.0, & 0.0, 0.0, 1.0, 45.0, &
0.0, 0.0, 1.0, 45.0 & 0.0, 0.0, 1.0, 45.0 &
],shape(CFTOCI_BAINROT)) ],shape(CFTOCI_BAINROT)),pReal)
if (present(a_cI) .and. present(a_cF)) then if (present(a_cI) .and. present(a_cF)) then
do i = 1,sum(Ntrans) do i = 1,sum(Ntrans)

4
src/math.f90
View File

@ -1110,7 +1110,7 @@ pure function math_rotationalPart(F) result(R)
C = matmul(transpose(F),F) C = matmul(transpose(F),F)
I_C = math_invariantsSym33(C) I_C = math_invariantsSym33(C)
I_F = [math_trace33(F), 0.5*(math_trace33(F)**2 - math_trace33(matmul(F,F)))] I_F = [math_trace33(F), 0.5_pReal*(math_trace33(F)**2 - math_trace33(matmul(F,F)))]
x = math_clip(I_C(1)**2 -3.0_pReal*I_C(2),0.0_pReal)**(3.0_pReal/2.0_pReal) x = math_clip(I_C(1)**2 -3.0_pReal*I_C(2),0.0_pReal)**(3.0_pReal/2.0_pReal)
if (dNeq0(x)) then if (dNeq0(x)) then
@ -1386,7 +1386,7 @@ subroutine selfTest()
call random_number(v3_3) call random_number(v3_3)
call random_number(v3_4) call random_number(v3_4)
if (dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, & if (dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0_pReal, &
math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) & math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) &
error stop 'math_volTetrahedron' error stop 'math_volTetrahedron'

2
src/mesh/FEM_utilities.f90
View File

@ -131,7 +131,7 @@ subroutine FEM_utilities_init
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),err_PETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal) wgt = real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)**(-1)
end subroutine FEM_utilities_init end subroutine FEM_utilities_init

2
src/mesh/discretization_mesh.f90
View File

@ -247,7 +247,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI) mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI)
do dirJ = 1, dimPlex do dirJ = 1, dimPlex
mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + & mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + &
pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0) pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0_pReal)
enddo enddo
enddo enddo
qOffset = qOffset + dimPlex qOffset = qOffset + dimPlex

102
src/mesh/mesh_mech_FEM.f90
View File

@ -230,14 +230,14 @@ subroutine FEM_mechanical_init(fieldBC)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetConvergenceTest(mechanical_snes,FEM_mechanical_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,err_PETSc) call SNESSetConvergenceTest(mechanical_snes,FEM_mechanical_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetTolerances(mechanical_snes,1.0,0.0,0.0,num%itmax,num%itmax,err_PETSc) call SNESSetTolerances(mechanical_snes,1.0_pReal,0.0_pReal,0.0_pReal,num%itmax,num%itmax,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetFromOptions(mechanical_snes,err_PETSc); CHKERRQ(err_PETSc) call SNESSetFromOptions(mechanical_snes,err_PETSc); CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call VecSet(solution ,0.0,err_PETSc); CHKERRQ(err_PETSc) call VecSet(solution ,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
call VecSet(solution_rate ,0.0,err_PETSc); CHKERRQ(err_PETSc) call VecSet(solution_rate,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
allocate(x_scal(cellDof)) allocate(x_scal(cellDof))
allocate(nodalWeightsP(1)) allocate(nodalWeightsP(1))
allocate(nodalPointsP(dimPlex)) allocate(nodalPointsP(dimPlex))
@ -338,11 +338,10 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx, & qPt, basis, comp, cidx, &
numFields, & numFields, &
bcSize,m bcSize,m,i
PetscReal :: detFAvg, detJ PetscReal :: detFAvg, detJ
PetscReal, dimension(dimPlex*dimPlex,cellDof) :: BMat PetscReal, dimension(dimPlex*dimPlex,cellDof) :: BMat
IS :: bcPoints
IS :: bcPoints
allocate(pV0(dimPlex)) allocate(pV0(dimPlex))
@ -358,7 +357,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMGetLocalVector(dm_local,x_local,err_PETSc) call DMGetLocalVector(dm_local,x_local,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecWAXPY(x_local,1.0,xx_local,solution_local,err_PETSc) call VecWAXPY(x_local,1.0_pReal,xx_local,solution_local,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then if (params%fieldBC%componentBC(field)%Mask(face)) then
@ -386,23 +385,23 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
do qPt = 0, nQuadrature-1 do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1 m = cell*nQuadrature + qPt+1
BMat = 0.0 BMat = 0.0_pReal
do basis = 0, nBasis-1 do basis = 0, nBasis-1
do comp = 0, dimPlex-1 do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = & i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: & BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) matmul(IcellJMat,basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex))
enddo enddo
enddo enddo
homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
enddo enddo
if (num%BBarStabilisation) then if (num%BBarStabilisation) then
detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature)) detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature,pReal))
do qPt = 0, nQuadrature-1 do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1 m = cell*nQuadrature + qPt+1
homogenization_F(1:dimPlex,1:dimPlex,m) = homogenization_F(1:dimPlex,1:dimPlex,m) & homogenization_F(1:dimPlex,1:dimPlex,m) = homogenization_F(1:dimPlex,1:dimPlex,m) &
* (detFAvg/math_det33(homogenization_F(1:3,1:3,m)))**(1.0/real(dimPlex)) * (detFAvg/math_det33(homogenization_F(1:3,1:3,m)))**(1.0_pReal/real(dimPlex,pReal))
enddo enddo
endif endif
@ -425,22 +424,22 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc) call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
f_scal = 0.0 f_scal = 0.0_pReal
do qPt = 0, nQuadrature-1 do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1 m = cell*nQuadrature + qPt+1
BMat = 0.0 BMat = 0.0_pReal
do basis = 0, nBasis-1 do basis = 0, nBasis-1
do comp = 0, dimPlex-1 do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = & i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: & BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) matmul(IcellJMat,basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex))
enddo enddo
enddo enddo
f_scal = f_scal + & f_scal = f_scal &
matmul(transpose(BMat), & + matmul(transpose(BMat), &
reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,m)), & reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,m)), &
shape=[dimPlex*dimPlex]))*qWeights(qPt+1) shape=[dimPlex*dimPlex]))*qWeights(qPt+1_pPETSCINT)
enddo enddo
f_scal = f_scal*abs(detJ) f_scal = f_scal*abs(detJ)
pf_scal => f_scal pf_scal => f_scal
@ -465,28 +464,25 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
PetscObject, intent(in) :: dummy PetscObject, intent(in) :: dummy
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
PetscDS :: prob PetscDS :: prob
Vec :: x_local, xx_local Vec :: x_local, xx_local
PetscSection :: section, gSection
PetscSection :: section, gSection
PetscReal, dimension(1, cellDof) :: MatB PetscReal, dimension(1, cellDof) :: MatB
PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA
PetscReal, dimension(3,3) :: F, FAvg, FInv PetscReal, dimension(3,3) :: F, FAvg, FInv
PetscReal :: detJ PetscReal :: detJ
PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & PetscReal, dimension(:), pointer :: basisField, basisFieldDer, &
pV0, pCellJ, pInvcellJ pV0, pCellJ, pInvcellJ
PetscScalar, dimension(:), pointer :: pK_e, x_scal PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscScalar,dimension(cellDOF,cellDOF), target :: K_e PetscScalar,dimension(cellDOF,cellDOF), target :: K_e
PetscScalar,dimension(cellDOF,cellDOF) :: K_eA , & PetscScalar,dimension(cellDOF,cellDOF) :: K_eA, K_eB
K_eB
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx,bcSize, m qPt, basis, comp, cidx,bcSize, m, i
IS :: bcPoints
IS :: bcPoints
allocate(pV0(dimPlex)) allocate(pV0(dimPlex))
@ -530,30 +526,29 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc) call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
K_eA = 0.0 K_eA = 0.0_pReal
K_eB = 0.0 K_eB = 0.0_pReal
MatB = 0.0 MatB = 0.0_pReal
FAvg = 0.0 FAvg = 0.0_pReal
BMatAvg = 0.0 BMatAvg = 0.0_pReal
do qPt = 0, nQuadrature-1 do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt + 1 m = cell*nQuadrature + qPt + 1
BMat = 0.0 BMat = 0.0_pReal
do basis = 0, nBasis-1 do basis = 0, nBasis-1
do comp = 0, dimPlex-1 do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = & i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
matmul( reshape(pInvcellJ, shape = [dimPlex,dimPlex]),& BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: & matmul(reshape(pInvcellJ,[dimPlex,dimPlex]),basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex))
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo enddo
enddo enddo
MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), & MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), &
shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1) shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1_pPETSCINT)
if (num%BBarStabilisation) then if (num%BBarStabilisation) then
F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex]) F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex])
FInv = math_inv33(F) FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex)) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0_pReal/real(dimPlex,pReal))
K_eB = K_eB - & K_eB = K_eB - &
matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex**2,1_pPETSCINT]), & matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex**2,1_pPETSCINT]), &
matmul(reshape(FInv(1:dimPlex,1:dimPlex), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
@ -568,10 +563,10 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
enddo enddo
if (num%BBarStabilisation) then if (num%BBarStabilisation) then
FInv = math_inv33(FAvg) FInv = math_inv33(FAvg)
K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + & K_e = K_eA*math_det33(FAvg/real(nQuadrature,pReal))**(1.0_pReal/real(dimPlex,pReal)) + &
(matmul(matmul(transpose(BMatAvg), & (matmul(matmul(transpose(BMatAvg), &
reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex**2,1_pPETSCINT],order=[2,1])),MatB) + & reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex**2,1_pPETSCINT],order=[2,1])),MatB) + &
K_eB)/real(dimPlex) K_eB)/real(dimPlex,pReal)
else else
K_e = K_eA K_e = K_eA
endif endif
@ -641,7 +636,7 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,err_PETSc) call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecAXPY(solution_local,1.0,x_local,err_PETSc); CHKERRQ(err_PETSc) call VecAXPY(solution_local,1.0_pReal,x_local,err_PETSc); CHKERRQ(err_PETSc)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(face)) then if (fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
@ -659,7 +654,7 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update rate and forward last inc ! update rate and forward last inc
call VecCopy(solution,solution_rate,err_PETSc); CHKERRQ(err_PETSc) call VecCopy(solution,solution_rate,err_PETSc); CHKERRQ(err_PETSc)
call VecScale(solution_rate,1.0/timeinc_old,err_PETSc); CHKERRQ(err_PETSc) call VecScale(solution_rate,timeinc_old**(-1),err_PETSc); CHKERRQ(err_PETSc)
endif endif
call VecCopy(solution_rate,solution,err_PETSc); CHKERRQ(err_PETSc) call VecCopy(solution_rate,solution,err_PETSc); CHKERRQ(err_PETSc)
call VecScale(solution,timeinc,err_PETSc); CHKERRQ(err_PETSc) call VecScale(solution,timeinc,err_PETSc); CHKERRQ(err_PETSc)
@ -685,9 +680,8 @@ subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reaso
call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,err_PETSc) call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN
print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), & print'(/,1x,a,a,i0,a,f0.3)', trim(incInfo), &
' @ Iteration ',PETScIter,' mechanical residual norm = ', & ' @ Iteration ',PETScIter,' mechanical residual norm = ',fnorm/divTol
int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol) ! ToDo: int casting?
print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
'Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal 'Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal
flush(IO_STDOUT) flush(IO_STDOUT)

6
src/phase_mechanical.f90
View File

@ -593,7 +593,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b
iteration: do NiterationState = 1, num%nState iteration: do NiterationState = 1, num%nState
dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1) dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0_pReal, nIterationState > 1)
dotState_last(1:sizeDotState,1) = dotState dotState_last(1:sizeDotState,1) = dotState
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en) broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
@ -756,7 +756,7 @@ function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b
real(pReal), dimension(3), parameter :: & real(pReal), dimension(3), parameter :: &
C = [0.5_pReal, 0.5_pReal, 1.0_pReal] C = [0.5_pReal, 0.5_pReal, 1.0_pReal]
real(pReal), dimension(4), parameter :: & real(pReal), dimension(4), parameter :: &
B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] B = [6.0_pReal, 3.0_pReal, 3.0_pReal, 6.0_pReal]**(-1)
broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C) broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C)
@ -1263,8 +1263,6 @@ module subroutine mechanical_restartRead(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph integer, intent(in) :: ph
integer :: en
call HDF5_read(plasticState(ph)%state0,groupHandle,'omega_plastic') call HDF5_read(plasticState(ph)%state0,groupHandle,'omega_plastic')
call HDF5_read(phase_mechanical_S0(ph)%data,groupHandle,'S') call HDF5_read(phase_mechanical_S0(ph)%data,groupHandle,'S')

4
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -912,7 +912,7 @@ pure subroutine kinetics_tw(Mp,T,abs_dot_gamma_sl,ph,en,&
real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: &
ddot_gamma_dtau_tw ddot_gamma_dtau_tw
real :: & real(pReal) :: &
tau, tau_r, tau_hat, & tau, tau_r, tau_hat, &
dot_N_0, & dot_N_0, &
x0, V, & x0, V, &
@ -988,7 +988,7 @@ pure subroutine kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,&
real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: &
ddot_gamma_dtau_tr ddot_gamma_dtau_tr
real :: & real(pReal) :: &
tau, tau_r, tau_hat, & tau, tau_r, tau_hat, &
dot_N_0, & dot_N_0, &
x0, V, & x0, V, &

2
src/rotations.f90
View File

@ -801,7 +801,7 @@ subroutine selfTest()
real(pReal), dimension(3,3) :: om, t33 real(pReal), dimension(3,3) :: om, t33
real(pReal), dimension(3,3,3,3) :: t3333 real(pReal), dimension(3,3,3,3) :: t3333
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
real :: A,B real(pReal) :: A,B
integer :: i integer :: i