diff --git a/CMakeLists.txt b/CMakeLists.txt index d24d9b8b8..6827d31b8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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}") 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}") if(CMAKE_BUILD_TYPE STREQUAL "DEBUG") diff --git a/PRIVATE b/PRIVATE index b14f78e96..a561f30a9 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit b14f78e96a8e2986aaf6845b98ea77fec92bc997 +Subproject commit a561f30a96c5b90e220b905a24c9d54acd479497 diff --git a/VERSION b/VERSION index 12c5e4193..45349346b 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha6-476-gac626db9d +3.0.0-alpha6-483-gde89d0fe3 diff --git a/cmake/Compiler-GNU.cmake b/cmake/Compiler-GNU.cmake index 6eedffc21..397d1277b 100644 --- a/cmake/Compiler-GNU.cmake +++ b/cmake/Compiler-GNU.cmake @@ -135,10 +135,3 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fsanitize=undefined") # detect undefined behavior # Additional options # -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 diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 7f34e4a13..4125aa8ef 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -118,8 +118,3 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all") # -check: Checks at runtime, where # arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays (performance?) # 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) diff --git a/cmake/Compiler-IntelLLVM.cmake b/cmake/Compiler-IntelLLVM.cmake index 883873e1c..bd0f07ee8 100644 --- a/cmake/Compiler-IntelLLVM.cmake +++ b/cmake/Compiler-IntelLLVM.cmake @@ -117,8 +117,3 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all") # -check: Checks at runtime, where # arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays (performance?) # 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) diff --git a/src/Marc/DAMASK_Marc.f90 b/src/Marc/DAMASK_Marc.f90 index 01eb8c5ef..c1f5b1d68 100644 --- a/src/Marc/DAMASK_Marc.f90 +++ b/src/Marc/DAMASK_Marc.f90 @@ -211,7 +211,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & use OMP_LIB implicit none(type,external) - 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 nn, & !< integration point number 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 ifr, & !< set to 1 if R 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 matus, & !< (1) user material identification number, (2) internal material identification 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 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) s = stress(1:ndi+nshear) @@ -383,16 +383,17 @@ subroutine flux(f,ts,n,time) use discretization_Marc implicit none(type,external) - real(pReal), dimension(6), intent(in) :: & + real(pReal), dimension(6), intent(in) :: & ts - integer, dimension(10), intent(in) :: & + integer(pI64), dimension(10), intent(in) :: & n - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & time - real(pReal), dimension(2), intent(out) :: & + real(pReal), dimension(2), intent(out) :: & 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 end subroutine flux @@ -410,7 +411,8 @@ subroutine uedinc(inc,incsub) use discretization_Marc implicit none(type,external) - integer, intent(in) :: inc, incsub + integer(pI64), intent(in) :: inc, incsub + integer :: n, nqncomp, nqdatatype integer, save :: inc_written real(pReal), allocatable, dimension(:,:) :: d_n @@ -427,9 +429,9 @@ subroutine uedinc(inc,incsub) enddo 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 end subroutine uedinc diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index dbaa668c2..4395c1581 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -276,7 +276,7 @@ program DAMASK_grid write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'R:',& 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)%N < 1) errorID = 835 if (loadCases(l)%f_out < 1) errorID = 836 @@ -505,7 +505,7 @@ subroutine getMaskedTensor(values,mask,tensor) integer :: i,j - values = 0.0 + values = 0.0_pReal do i = 1,3 row => tensor%get(i) do j = 1,3 diff --git a/src/grid/VTI.f90 b/src/grid/VTI.f90 index 98cb55e54..72f57fb86 100644 --- a/src/grid/VTI.f90 +++ b/src/grid/VTI.f90 @@ -151,7 +151,7 @@ subroutine VTI_readCellsSizeOrigin(cells,geomSize,origin, & character(len=*), intent(in) :: & fileContent - character(len=:), allocatable :: dataType, headerType + character(len=:), allocatable :: headerType logical :: inFile, inImage, compressed integer(pI64) :: & startPos, endPos diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 682a3e5b2..158ee0a8d 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -233,9 +233,9 @@ pure function cellSurfaceArea(geomSize,cells) real(pReal), dimension(6,1,product(cells)) :: cellSurfaceArea - cellSurfaceArea(1:2,1,:) = geomSize(2)/real(cells(2)) * geomSize(3)/real(cells(3)) - cellSurfaceArea(3:4,1,:) = geomSize(3)/real(cells(3)) * geomSize(1)/real(cells(1)) - cellSurfaceArea(5:6,1,:) = geomSize(1)/real(cells(1)) * geomSize(2)/real(cells(2)) + 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),pReal) * geomSize(1)/real(cells(1),pReal) + cellSurfaceArea(5:6,1,:) = geomSize(1)/real(cells(1),pReal) * geomSize(2)/real(cells(2),pReal) end function cellSurfaceArea diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 1040d6fb7..4bb705807 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -226,14 +226,14 @@ subroutine grid_mechanical_FEM_init delta = geomSize/real(cells,pReal) ! grid spacing detJ = product(delta) ! cell volume - BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & - 1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & - -1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & - 1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & - -1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & - 1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & - -1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3), & - 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 + BMat = reshape(real([-delta(1)**(-1),-delta(2)**(-1),-delta(3)**(-1), & + delta(1)**(-1),-delta(2)**(-1),-delta(3)**(-1), & + -delta(1)**(-1), delta(2)**(-1),-delta(3)**(-1), & + delta(1)**(-1), delta(2)**(-1),-delta(3)**(-1), & + -delta(1)**(-1),-delta(2)**(-1), delta(3)**(-1), & + delta(1)**(-1),-delta(2)**(-1), delta(3)**(-1), & + -delta(1)**(-1), delta(2)**(-1), delta(3)**(-1), & + delta(1)**(-1), delta(2)**(-1), delta(3)**(-1)],pReal), [3,8])/4.0_pReal ! shape function derivative matrix HGMat = matmul(transpose(HGcomp),HGcomp) & * HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal ! hourglass stabilization matrix @@ -661,7 +661,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc) MatNullSpace :: matnull PetscErrorCode :: err_PETSc - BMatFull = 0.0 + BMatFull = 0.0_pReal BMatFull(1:3,1 :8 ) = BMat BMatFull(4:6,9 :16) = BMat BMatFull(7:9,17:24) = BMat @@ -691,7 +691,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc) enddo; enddo; enddo row = col 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) + & homogenization_dPdF(2,2,2,2,ce) + & homogenization_dPdF(3,3,3,3,ce))/3.0_pReal diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index aa2ed1de9..7b6225577 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -203,7 +203,7 @@ subroutine spectral_utilities_init() CHKERRQ(err_PETSc) 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%divergence_correction = num_grid%get_asInt('divergence_correction', defaultVal=2) @@ -270,7 +270,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],& 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) 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]) @@ -279,7 +279,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],& 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) 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]) @@ -288,7 +288,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,& 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) call c_f_pointer(scalarField,scalarField_real, & [2_C_INTPTR_T*(cellsFFTW(1)/2_C_INTPTR_T+1_C_INTPTR_T),cellsFFTW(2),z]) @@ -390,17 +390,17 @@ subroutine utilities_updateGamma(C) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset) end do do concurrent(l = 1:3, m = 1:3) - temp33_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 #else forall(l = 1:3, m = 1:3) & xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset) 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 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 - 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) temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal) #ifndef __INTEL_COMPILER @@ -523,17 +523,17 @@ subroutine utilities_fourierGammaConvolution(fieldAim) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k) end do 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 #else 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_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 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 - 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) temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal) #ifndef __INTEL_COMPILER @@ -592,8 +592,8 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t) !$OMP PARALLEL DO PRIVATE(GreenOp_hat) do k = 1, cells3; do j = 1, cells(2) ;do i = 1, cells1Red GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) & - / (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal) & - * sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k)))) + / (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,pReal),xi1st(1:3,i,j,k)))) scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat enddo; enddo; enddo !$OMP END PARALLEL DO @@ -613,7 +613,7 @@ real(pReal) function utilities_divergenceRMS() print'(/,1x,a)', '... calculating divergence ................................................' 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 @@ -657,7 +657,7 @@ real(pReal) function utilities_curlRMS() print'(/,1x,a)', '... calculating curl ......................................................' 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 diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 2e66da591..3363408d0 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -3,8 +3,6 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization) damage - use lattice - interface module subroutine pass_init @@ -79,8 +77,9 @@ end subroutine damage_init !-------------------------------------------------------------------------------------------------- module subroutine damage_partition(ce) + integer, intent(in) :: ce + real(pReal) :: phi - integer, intent(in) :: ce if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return @@ -91,7 +90,7 @@ end subroutine damage_partition !-------------------------------------------------------------------------------------------------- -!> @brief Homogenized damage viscosity. +!> @brief Homogenize damage viscosity. !-------------------------------------------------------------------------------------------------- module function homogenization_mu_phi(ce) result(mu) @@ -105,7 +104,7 @@ end function homogenization_mu_phi !-------------------------------------------------------------------------------------------------- -!> @brief Homogenized damage conductivity/diffusivity in reference configuration. +!> @brief Homogenize damage conductivity. !-------------------------------------------------------------------------------------------------- module function homogenization_K_phi(ce) result(K) @@ -119,13 +118,12 @@ end function homogenization_K_phi !-------------------------------------------------------------------------------------------------- -!> @brief Homogenized damage driving force. +!> @brief Homogenize damage driving force. !-------------------------------------------------------------------------------------------------- module function homogenization_f_phi(phi,ce) result(f) integer, intent(in) :: ce - real(pReal), intent(in) :: & - phi + real(pReal), intent(in) :: phi real(pReal) :: f @@ -140,8 +138,7 @@ end function homogenization_f_phi module subroutine homogenization_set_phi(phi,ce) integer, intent(in) :: ce - real(pReal), intent(in) :: & - phi + real(pReal), intent(in) :: phi integer :: & ho, & diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index 911def035..8e8ae1df9 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -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 vPen = 0.0_pReal do i = 1,nGrain - vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* & - sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* & - gVol(i)*transpose(math_inv33(fDef(:,:,i))) + vPen(:,:,i) = -real(nGrain,pReal)**(-1)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr & + * sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0_pReal),vDiscrep) & + * gVol(i)*transpose(math_inv33(fDef(:,:,i))) end do end subroutine volumePenalty diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index e9c9d2195..ceed47365 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -3,8 +3,6 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization) thermal - use lattice - interface module subroutine pass_init @@ -89,7 +87,7 @@ end subroutine thermal_init !-------------------------------------------------------------------------------------------------- module subroutine thermal_partition(ce) - integer, intent(in) :: ce + integer, intent(in) :: ce real(pReal) :: T, dot_T integer :: co @@ -105,7 +103,7 @@ end subroutine thermal_partition !-------------------------------------------------------------------------------------------------- -!> @brief Homogenized thermal viscosity. +!> @brief Homogenize thermal viscosity. !-------------------------------------------------------------------------------------------------- module function homogenization_mu_T(ce) result(mu) @@ -124,7 +122,7 @@ end function homogenization_mu_T !-------------------------------------------------------------------------------------------------- -!> @brief Homogenized thermal conductivity in reference configuration. +!> @brief Homogenize thermal conductivity. !-------------------------------------------------------------------------------------------------- module function homogenization_K_T(ce) result(K) @@ -143,7 +141,7 @@ end function homogenization_K_T !-------------------------------------------------------------------------------------------------- -!> @brief Homogenized heat generation rate. +!> @brief Homogenize heat generation rate. !-------------------------------------------------------------------------------------------------- module function homogenization_f_T(ce) result(f) @@ -167,7 +165,7 @@ end function homogenization_f_T module subroutine homogenization_thermal_setField(T,dot_T, ce) integer, intent(in) :: ce - real(pReal), intent(in) :: T, dot_T + real(pReal), intent(in) :: T, dot_T current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce)) = T @@ -187,6 +185,7 @@ module subroutine thermal_results(ho,group) integer :: o + associate(prm => param(ho)) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) diff --git a/src/lattice.f90 b/src/lattice.f90 index 5339dc9ac..7bd846d0f 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2008,7 +2008,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI) ],pReal),shape(CFTOHP_SYSTEMTRANS)) 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, & 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, & 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 :: & CFTOCI_BAINVARIANT = reshape( [& @@ -2040,7 +2040,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI) ],shape(CFTOCI_BAINVARIANT)) 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, & 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 & - ],shape(CFTOCI_BAINROT)) + ],shape(CFTOCI_BAINROT)),pReal) if (present(a_cI) .and. present(a_cF)) then do i = 1,sum(Ntrans) diff --git a/src/math.f90 b/src/math.f90 index d11ad1565..07b657cab 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1110,7 +1110,7 @@ pure function math_rotationalPart(F) result(R) C = matmul(transpose(F),F) 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) if (dNeq0(x)) then @@ -1386,7 +1386,7 @@ subroutine selfTest() call random_number(v3_3) 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)) & error stop 'math_volTetrahedron' diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 4d24a10b3..3bdf901e5 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -140,7 +140,7 @@ subroutine FEM_utilities_init call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),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 diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index ad358e40b..bbb125887 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -253,7 +253,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints) mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI) do dirJ = 1, dimPlex 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 qOffset = qOffset + dimPlex diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 90c3aa31e..fae906891 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -243,14 +243,14 @@ subroutine FEM_mechanical_init(fieldBC) CHKERRQ(err_PETSc) call SNESSetConvergenceTest(mechanical_snes,FEM_mechanical_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,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) call SNESSetFromOptions(mechanical_snes,err_PETSc); CHKERRQ(err_PETSc) !-------------------------------------------------------------------------------------------------- ! init fields - call VecSet(solution ,0.0,err_PETSc); CHKERRQ(err_PETSc) - call VecSet(solution_rate ,0.0,err_PETSc); CHKERRQ(err_PETSc) + call VecSet(solution ,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc) + call VecSet(solution_rate,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc) allocate(x_scal(cellDof)) allocate(nodalWeightsP(1)) allocate(nodalPointsP(dimPlex)) @@ -351,11 +351,10 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc PetscInt :: cellStart, cellEnd, cell, field, face, & qPt, basis, comp, cidx, & numFields, & - bcSize,m + bcSize,m,i PetscReal :: detFAvg, detJ PetscReal, dimension(dimPlex*dimPlex,cellDof) :: BMat - - IS :: bcPoints + IS :: bcPoints allocate(pV0(dimPlex)) @@ -371,7 +370,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc CHKERRQ(err_PETSc) call DMGetLocalVector(dm_local,x_local,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) do field = 1, dimPlex; do face = 1, mesh_Nboundaries if (params%fieldBC%componentBC(field)%Mask(face)) then @@ -399,23 +398,23 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) do qPt = 0, nQuadrature-1 m = cell*nQuadrature + qPt+1 - BMat = 0.0 + BMat = 0.0_pReal do basis = 0, nBasis-1 do comp = 0, dimPlex-1 cidx = basis*dimPlex+comp - BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = & - matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: & - (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) + i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp + BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = & + matmul(IcellJMat,basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex)) enddo enddo homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) enddo 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 m = cell*nQuadrature + qPt+1 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 endif @@ -438,22 +437,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) CHKERRQ(err_PETSc) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) - f_scal = 0.0 + f_scal = 0.0_pReal do qPt = 0, nQuadrature-1 m = cell*nQuadrature + qPt+1 - BMat = 0.0 + BMat = 0.0_pReal do basis = 0, nBasis-1 do comp = 0, dimPlex-1 cidx = basis*dimPlex+comp - BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = & - matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: & - (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) + i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp + BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = & + matmul(IcellJMat,basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex)) enddo enddo - f_scal = f_scal + & - matmul(transpose(BMat), & - reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,m)), & - shape=[dimPlex*dimPlex]))*qWeights(qPt+1) + f_scal = f_scal & + + matmul(transpose(BMat), & + reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,m)), & + shape=[dimPlex*dimPlex]))*qWeights(qPt+1_pPETSCINT) enddo f_scal = f_scal*abs(detJ) pf_scal => f_scal @@ -478,28 +477,25 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P PetscObject, intent(in) :: dummy PetscErrorCode :: err_PETSc - PetscDS :: prob - Vec :: x_local, xx_local - - PetscSection :: section, gSection + PetscDS :: prob + Vec :: x_local, xx_local + PetscSection :: section, gSection PetscReal, dimension(1, cellDof) :: MatB PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA - PetscReal, dimension(3,3) :: F, FAvg, FInv - PetscReal :: detJ - PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & + PetscReal, dimension(3,3) :: F, FAvg, FInv + PetscReal :: detJ + PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & pV0, pCellJ, pInvcellJ PetscScalar, dimension(:), pointer :: pK_e, x_scal PetscScalar,dimension(cellDOF,cellDOF), target :: K_e - PetscScalar,dimension(cellDOF,cellDOF) :: K_eA , & - K_eB + PetscScalar,dimension(cellDOF,cellDOF) :: K_eA, K_eB - PetscInt :: cellStart, cellEnd, cell, field, face, & - qPt, basis, comp, cidx,bcSize, m - - IS :: bcPoints + PetscInt :: cellStart, cellEnd, cell, field, face, & + qPt, basis, comp, cidx,bcSize, m, i + IS :: bcPoints allocate(pV0(dimPlex)) @@ -543,30 +539,29 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P CHKERRQ(err_PETSc) call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc) CHKERRQ(err_PETSc) - K_eA = 0.0 - K_eB = 0.0 - MatB = 0.0 - FAvg = 0.0 - BMatAvg = 0.0 + K_eA = 0.0_pReal + K_eB = 0.0_pReal + MatB = 0.0_pReal + FAvg = 0.0_pReal + BMatAvg = 0.0_pReal do qPt = 0, nQuadrature-1 m = cell*nQuadrature + qPt + 1 - BMat = 0.0 + BMat = 0.0_pReal do basis = 0, nBasis-1 do comp = 0, dimPlex-1 cidx = basis*dimPlex+comp - BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = & - matmul( reshape(pInvcellJ, shape = [dimPlex,dimPlex]),& - basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: & - (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) + i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp + BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = & + matmul(reshape(pInvcellJ,[dimPlex,dimPlex]),basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex)) enddo enddo 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]),BMat)*qWeights(qPt+1) + shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1_pPETSCINT) if (num%BBarStabilisation) then F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex]) 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 - & matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex**2,1_pPETSCINT]), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), & @@ -581,10 +576,10 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P enddo if (num%BBarStabilisation) then 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), & 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 K_e = K_eA endif @@ -654,7 +649,7 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC) CHKERRQ(err_PETSc) call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,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 if (fieldBC%componentBC(field)%Mask(face)) then call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) @@ -672,7 +667,7 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC) !-------------------------------------------------------------------------------------------------- ! update rate and forward last inc 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 call VecCopy(solution_rate,solution,err_PETSc); CHKERRQ(err_PETSc) call VecScale(solution,timeinc,err_PETSc); CHKERRQ(err_PETSc) @@ -698,9 +693,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) CHKERRQ(err_PETSc) if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN - print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), & - ' @ Iteration ',PETScIter,' mechanical residual norm = ', & - int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol) ! ToDo: int casting? + print'(/,1x,a,a,i0,a,f0.3)', trim(incInfo), & + ' @ Iteration ',PETScIter,' mechanical residual norm = ',fnorm/divTol 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 flush(IO_STDOUT) diff --git a/src/phase.f90 b/src/phase.f90 index f9f3a8125..f304707b9 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -626,9 +626,10 @@ function crystallite_push33ToRef(co,ce, tensor33) ce real(pReal), dimension(3,3) :: crystallite_push33ToRef - real(pReal), dimension(3,3) :: T + real(pReal), dimension(3,3) :: T integer :: ph, en + ph = material_phaseID(co,ce) en = material_phaseEntry(co,ce) T = matmul(phase_O_0(ph)%data(en)%asMatrix(),transpose(math_inv33(phase_F(co,ce)))) ! ToDo: initial orientation correct? diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 82d94150d..7211f2db5 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -4,8 +4,9 @@ submodule(phase) damage type :: tDamageParameters - real(pReal) :: mu = 0.0_pReal !< viscosity - real(pReal), dimension(3,3) :: D = 0.0_pReal !< conductivity/diffusivity + real(pReal) :: & + mu = 0.0_pReal, & !< viscosity + l_c = 0.0_pReal !< characteristic length end type tDamageParameters enum, bind(c); enumerator :: & @@ -104,8 +105,8 @@ module subroutine damage_init if (sources%length == 1) then damage_active = .true. source => sources%get(1) - param(ph)%mu = source%get_asFloat('mu') - param(ph)%D = math_I3 * source%get_asFloat('l_c')**2 + param(ph)%mu = source%get_asFloat('mu') + param(ph)%l_c = source%get_asFloat('l_c') end if end do @@ -117,7 +118,7 @@ module subroutine damage_init where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID end if - phase_damage_maxSizeDotState = maxval(damageState%sizeDotState) + phase_damage_maxSizeDotState = maxval(damageState%sizeDotState) end subroutine damage_init @@ -157,9 +158,9 @@ module function phase_damage_C66(C66,ph,en) result(C66_degraded) damageType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) damageType - C66_degraded = C66 * damage_phi(ph,en)**2 + C66_degraded = C66 * damage_phi(ph,en)**2 case default damageType - C66_degraded = C66 + C66_degraded = C66 end select damageType end function phase_damage_C66 @@ -385,7 +386,7 @@ module function phase_K_phi(co,ce) result(K) real(pReal), dimension(3,3) :: K - K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%D) + K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%l_c**2*math_I3) end function phase_K_phi diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 93357bee5..3a66bdf95 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -593,7 +593,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b 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 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 :: & C = [0.5_pReal, 0.5_pReal, 1.0_pReal] 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) @@ -1263,8 +1263,6 @@ module subroutine mechanical_restartRead(groupHandle,ph) integer(HID_T), intent(in) :: groupHandle integer, intent(in) :: ph - integer :: en - call HDF5_read(plasticState(ph)%state0,groupHandle,'omega_plastic') call HDF5_read(phase_mechanical_S0(ph)%data,groupHandle,'S') diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index dfd84ffc9..f01a7e95d 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -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) :: & ddot_gamma_dtau_tw - real :: & + real(pReal) :: & tau, tau_r, tau_hat, & dot_N_0, & 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) :: & ddot_gamma_dtau_tr - real :: & + real(pReal) :: & tau, tau_r, tau_hat, & dot_N_0, & x0, V, & diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 0f5836df9..ed374e142 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -4,8 +4,8 @@ submodule(phase) thermal type :: tThermalParameters - real(pReal) :: C_p = 0.0_pReal !< heat capacity - real(pReal), dimension(3,3) :: K = 0.0_pReal !< thermal conductivity + real(pReal) :: C_p = 0.0_pReal !< heat capacity + real(pReal), dimension(3,3) :: K = 0.0_pReal !< thermal conductivity character(len=pStringLen), allocatable, dimension(:) :: output end type tThermalParameters @@ -72,7 +72,7 @@ submodule(phase) thermal contains !---------------------------------------------------------------------------------------------- -!< @brief initializes thermal sources and kinematics mechanism +!< @brief Initializes thermal sources and kinematics mechanism. !---------------------------------------------------------------------------------------------- module subroutine thermal_init(phases) @@ -146,7 +146,7 @@ end subroutine thermal_init !---------------------------------------------------------------------------------------------- -!< @brief calculates thermal dissipation rate +!< @brief Calculate thermal source. !---------------------------------------------------------------------------------------------- module function phase_f_T(ph,en) result(f) @@ -176,7 +176,7 @@ end function phase_f_T !-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!> @brief tbd. !-------------------------------------------------------------------------------------------------- function phase_thermal_collectDotState(ph,en) result(broken) @@ -216,7 +216,7 @@ end function phase_mu_T !-------------------------------------------------------------------------------------------------- -!> @brief Thermal conductivity/diffusivity in reference configuration. +!> @brief Thermal conductivity in reference configuration. !-------------------------------------------------------------------------------------------------- module function phase_K_T(co,ce) result(K) @@ -255,6 +255,7 @@ function integrateThermalState(Delta_t, ph,en) result(broken) so, & sizeDotState + broken = phase_thermal_collectDotState(ph,en) if (broken) return diff --git a/src/rotations.f90 b/src/rotations.f90 index 83d945ea9..5c8677c81 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -801,7 +801,7 @@ subroutine selfTest() real(pReal), dimension(3,3) :: om, t33 real(pReal), dimension(3,3,3,3) :: t3333 real(pReal), dimension(6,6) :: C - real :: A,B + real(pReal) :: A,B integer :: i