diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a0cd44faa..b9b5fafff 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -176,6 +176,7 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral") add_library(SPECTRAL_SOLVER OBJECT "grid_thermal_spectral.f90" "grid_damage_spectral.f90" + "grid_mech_FEM.f90" "grid_mech_spectral_basic.f90" "grid_mech_spectral_polarisation.f90") add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES) diff --git a/src/grid_mech_FEM.f90 b/src/grid_mech_FEM.f90 new file mode 100644 index 000000000..029c05b4e --- /dev/null +++ b/src/grid_mech_FEM.f90 @@ -0,0 +1,735 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Arko Jyoti Bhattacharjee, Max-Planck-Institut für Eisenforschung GmbH +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Grid solver for mechanics: FEM +!-------------------------------------------------------------------------------------------------- +module grid_mech_FEM +#include +#include + use PETScdmda + use PETScsnes + use prec, only: & + pInt, & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private + + character (len=*), parameter, public :: & + grid_mech_FEM_label = 'fem' + +!-------------------------------------------------------------------------------------------------- +! derived types + type(tSolutionParams), private :: params + +!-------------------------------------------------------------------------------------------------- +! PETSc data + DM, private :: mech_grid + SNES, private :: mech_snes + Vec, private :: solution_current, solution_lastInc, solution_rate + +!-------------------------------------------------------------------------------------------------- +! common pointwise data + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc + real(pReal), private :: detJ + real(pReal), private, dimension(3) :: delta + real(pReal), private, dimension(3,8) :: BMat + real(pReal), private, dimension(8,8) :: HGMat + PetscInt, private :: xstart,ystart,zstart,xend,yend,zend + +!-------------------------------------------------------------------------------------------------- +! stress, stiffness and compliance average etc. + real(pReal), private, dimension(3,3) :: & + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient + F_aim = math_I3, & !< current prescribed deformation gradient + F_aim_lastIter = math_I3, & + F_aim_lastInc = math_I3, & !< previous average deformation gradient + P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress + + character(len=1024), private :: incInfo !< time and increment information + + real(pReal), private, dimension(3,3,3,3) :: & + C_volAvg = 0.0_pReal, & !< current volume average stiffness + C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness + S = 0.0_pReal !< current compliance (filled up with zeros) + + real(pReal), private :: & + err_BC !< deviation from stress BC + + integer(pInt), private :: & + totalIter = 0_pInt !< total iteration in current increment + + public :: & + grid_mech_FEM_init, & + grid_mech_FEM_solution, & + grid_mech_FEM_forward + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all necessary fields and fills them with data, potentially from restart info +!-------------------------------------------------------------------------------------------------- +subroutine grid_mech_FEM_init + use IO, only: & + IO_intOut, & + IO_error + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRestart + use FEsolving, only: & + restartInc + use numerics, only: & + worldrank, & + worldsize, & + petsc_options + use homogenization, only: & + materialpoint_F0 + use DAMASK_interface, only: & + getSolverJobName + use spectral_utilities, only: & + utilities_constitutiveResponse, & + utilities_updateIPcoords, & + wgt + use mesh, only: & + geomSize, & + grid, & + grid3 + use math, only: & + math_invSym3333 + + implicit none + real(pReal) :: HGCoeff = 0e-2_pReal + PetscInt, dimension(:), allocatable :: localK + real(pReal), dimension(3,3) :: & + temp33_Real = 0.0_pReal + real(pReal), dimension(4,8) :: & + HGcomp = reshape([ 1.0_pReal, 1.0_pReal, 1.0_pReal,-1.0_pReal, & + 1.0_pReal,-1.0_pReal,-1.0_pReal, 1.0_pReal, & + -1.0_pReal, 1.0_pReal,-1.0_pReal, 1.0_pReal, & + -1.0_pReal,-1.0_pReal, 1.0_pReal,-1.0_pReal, & + -1.0_pReal,-1.0_pReal, 1.0_pReal, 1.0_pReal, & + -1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, & + 1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, & + 1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8]) + PetscErrorCode :: ierr + integer(pInt) :: rank + character(len=1024) :: rankStr + real(pReal), dimension(3,3,3,3) :: devNull + PetscScalar, pointer, dimension(:,:,:,:) :: & + u_current,u_lastincrement,u_rate + write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>' + +!-------------------------------------------------------------------------------------------------- +! set default and user defined options for PETSc + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres & + &-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr) + CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! allocate global fields + allocate (F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (P_current (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + +!-------------------------------------------------------------------------------------------------- +! initialize solver specific parts of PETSc + call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) + allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 + do rank = 1, worldsize + call MPI_Bcast(localK(rank),1,MPI_INTEGER,rank-1,PETSC_COMM_WORLD,ierr) + enddo + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, & + DMDA_STENCIL_BOX, & + grid(1),grid(2),grid(3), & + 1, 1, worldsize, & + 3, 1, & + [grid(1)],[grid(2)],localK, & + mech_grid,ierr) + CHKERRQ(ierr) + call DMDASetUniformCoordinates(mech_grid,0.0,geomSize(1),0.0,geomSize(2),0.0,geomSize(3),ierr) + CHKERRQ(ierr) + call SNESSetDM(mech_snes,mech_grid,ierr); CHKERRQ(ierr) + call DMsetFromOptions(mech_grid,ierr); CHKERRQ(ierr) + call DMsetUp(mech_grid,ierr); CHKERRQ(ierr) + call DMCreateGlobalVector(mech_grid,solution_current,ierr); CHKERRQ(ierr) + call DMCreateGlobalVector(mech_grid,solution_lastInc,ierr); CHKERRQ(ierr) + call DMCreateGlobalVector(mech_grid,solution_rate ,ierr); CHKERRQ(ierr) + call DMSNESSetFunctionLocal(mech_grid,grid_mech_FEM_formResidual,PETSC_NULL_SNES,ierr) + CHKERRQ(ierr) + call DMSNESSetJacobianLocal(mech_grid,grid_mech_FEM_formJacobian,PETSC_NULL_SNES,ierr) + CHKERRQ(ierr) + call SNESSetConvergenceTest(mech_snes,grid_mech_FEM_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) + CHKERRQ(ierr) ! specify custom convergence check function "_converged" + call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures + call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments + +!-------------------------------------------------------------------------------------------------- +! init fields + call VecSet(solution_current,0.0,ierr);CHKERRQ(ierr) + call VecSet(solution_lastInc,0.0,ierr);CHKERRQ(ierr) + call VecSet(solution_rate ,0.0,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastincrement,ierr); CHKERRQ(ierr) + + call DMDAGetCorners(mech_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent + CHKERRQ(ierr) + xend = xstart+xend-1 + yend = ystart+yend-1 + zend = zstart+zend-1 + delta = geomSize/real(grid,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 + + HGMat = matmul(transpose(HGcomp),HGcomp) & + * HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal ! hourglass stabilization matrix + + restart: if (restartInc > 0_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading values of increment ', restartInc, ' from file' + flush(6) + endif + !write(rankStr,'(a1,i0)')'_',worldrank + !call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F)) + !read (777,rec=1) F; close (777) + !call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc)) + !read (777,rec=1) F_lastInc; close (777) + !call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot)) + !read (777,rec=1) F_aimDot; close (777) + !call IO_read_realFile(777,'u_current'//trim(rankStr),trim(getSolverJobName()),size(u_current)) + !read (777,rec=1) u_current; close (777) + !call IO_read_realFile(777,'u_lastincrement'//trim(rankStr),trim(getSolverJobName()),size(u_lastincrement)) + !read (777,rec=1) u_lastincrement; close (777) + + F_aim = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt ! average of F + call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim') + F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt! average of F_lastInc + call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim_lastInc') + elseif (restartInc == 0_pInt) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) + endif restart + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + call Utilities_updateIPcoords(F) + call Utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + F, & ! target F + 0.0_pReal, & ! time increment + math_I3) ! no rotation of boundary condition + call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) + CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastincrement,ierr) + CHKERRQ(ierr) + + restartRead: if (restartInc > 0_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0 .and. worldrank == 0_pInt) & + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading more values of increment ', restartInc, ' from file' + flush(6) + !call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) + !read (777,rec=1) C_volAvg; close (777) + !call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) + !read (777,rec=1) C_volAvgLastInc; close (777) + endif restartRead + +end subroutine grid_mech_FEM_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief solution for the FEM scheme with internal iterations +!-------------------------------------------------------------------------------------------------- +function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) + use IO, only: & + IO_error + use spectral_utilities, only: & + tBoundaryCondition, & + utilities_maskedCompliance + use FEsolving, only: & + restartWrite, & + terminallyIll + + implicit none + +!-------------------------------------------------------------------------------------------------- +! input data for solution + character(len=*), intent(in) :: & + incInfoIn + real(pReal), intent(in) :: & + timeinc, & !< time increment of current solution + timeinc_old !< time increment of last successful increment + type(tBoundaryCondition), intent(in) :: & + stress_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(tSolutionState) :: & + solution + +!-------------------------------------------------------------------------------------------------- +! PETSc Data + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + incInfo = incInfoIn + +!-------------------------------------------------------------------------------------------------- +! update stiffness (and gamma operator) + S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) +!-------------------------------------------------------------------------------------------------- +! set module wide available data + params%stress_mask = stress_BC%maskFloat + params%stress_BC = stress_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + +!-------------------------------------------------------------------------------------------------- +! solve BVP + call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr);CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! check convergence + call SNESGetConvergedReason(mech_snes,reason,ierr);CHKERRQ(ierr) + + solution%converged = reason > 0 + solution%iterationsNeeded = totalIter + solution%termIll = terminallyIll + terminallyIll = .false. + + if (reason == SNES_DIVERGED_FNORM_NAN) call IO_error(893_pInt) + +end function grid_mech_FEM_solution + + +!-------------------------------------------------------------------------------------------------- +!> @brief forms the residual vector +!-------------------------------------------------------------------------------------------------- +subroutine grid_mech_FEM_formResidual(da_local,x_local,f_local,dummy,ierr) + use numerics, only: & + itmax, & + itmin + use numerics, only: & + worldrank + use mesh, only: & + grid + use math, only: & + math_rotate_backward33, & + math_mul3333xx33 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation + use spectral_utilities, only: & + utilities_constitutiveResponse + use IO, only: & + IO_intOut + use FEsolving, only: & + terminallyIll + use homogenization, only: & + materialpoint_dPdF + + implicit none + DM :: da_local + Vec :: x_local, f_local + PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal + PetscScalar, dimension(8,3) :: x_elem, f_elem + PetscInt :: i, ii, j, jj, k, kk, ctr, ele + real(pReal), dimension(3,3) :: & + deltaF_aim + PetscInt :: & + PETScIter, & + nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal), dimension(3,3,3,3) :: devNull + + + call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) + + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment +!-------------------------------------------------------------------------------------------------- +! begin of new iteration + newIteration: if (totalIter <= PETScIter) then + totalIter = totalIter + 1_pInt + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & + trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim =', transpose(F_aim) + flush(6) + endif newIteration + +!-------------------------------------------------------------------------------------------------- +! get deformation gradient + call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) + enddo; enddo; enddo + ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 + F(1:3,1:3,ii,jj,kk) = F_aim + transpose(matmul(BMat,x_elem)) + enddo; enddo; enddo + call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! evaluate constitutive response + call Utilities_constitutiveResponse(P_current,& + P_av,C_volAvg,devNull, & + F,params%timeinc,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) + +!-------------------------------------------------------------------------------------------------- +! stress BC handling + F_aim_lastIter = F_aim + deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) + F_aim = F_aim - deltaF_aim + err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + +!-------------------------------------------------------------------------------------------------- +! constructing residual + call VecSet(f_local,0.0,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + ele = 0 + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) + enddo; enddo; enddo + ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 + ele = ele + 1 + f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + & + matmul(HGMat,x_elem)*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + f_scal(0:2,i+ii,j+jj,k+kk) = f_scal(0:2,i+ii,j+jj,k+kk) + f_elem(ctr,1:3) + enddo; enddo; enddo + enddo; enddo; enddo + call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! applying boundary conditions + call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + if (zstart == 0) then + f_scal(0:2,xstart,ystart,zstart) = 0.0 + f_scal(0:2,xend+1,ystart,zstart) = 0.0 + f_scal(0:2,xstart,yend+1,zstart) = 0.0 + f_scal(0:2,xend+1,yend+1,zstart) = 0.0 + endif + if (zend + 1 == grid(3)) then + f_scal(0:2,xstart,ystart,zend+1) = 0.0 + f_scal(0:2,xend+1,ystart,zend+1) = 0.0 + f_scal(0:2,xstart,yend+1,zend+1) = 0.0 + f_scal(0:2,xend+1,yend+1,zend+1) = 0.0 + endif + call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + +end subroutine grid_mech_FEM_formResidual + + +!-------------------------------------------------------------------------------------------------- +!> @brief forms the FEM stiffness matrix +!-------------------------------------------------------------------------------------------------- +subroutine grid_mech_FEM_formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr) + use mesh, only: & + mesh_ipCoordinates + use homogenization, only: & + materialpoint_dPdF + + implicit none + + DM :: da_local + Vec :: x_local, coordinates + Mat :: Jac_pre, Jac + MatStencil,dimension(4,24) :: row, col + PetscScalar,pointer,dimension(:,:,:,:) :: x_scal + PetscScalar,dimension(24,24) :: K_ele + PetscScalar,dimension(9,24) :: BMatFull + PetscInt :: i, ii, j, jj, k, kk, ctr, ele + PetscInt,dimension(3) :: rows + PetscScalar :: diag + PetscObject :: dummy + MatNullSpace :: matnull + PetscErrorCode :: ierr + + BMatFull = 0.0 + BMatFull(1:3,1 :8 ) = BMat + BMatFull(4:6,9 :16) = BMat + BMatFull(7:9,17:24) = BMat + call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr) + call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr) + call MatZeroEntries(Jac,ierr); CHKERRQ(ierr) + ele = 0 + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ctr = 0 + do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 + ctr = ctr + 1 + col(MatStencil_i,ctr ) = i+ii + col(MatStencil_j,ctr ) = j+jj + col(MatStencil_k,ctr ) = k+kk + col(MatStencil_c,ctr ) = 0 + col(MatStencil_i,ctr+8 ) = i+ii + col(MatStencil_j,ctr+8 ) = j+jj + col(MatStencil_k,ctr+8 ) = k+kk + col(MatStencil_c,ctr+8 ) = 1 + col(MatStencil_i,ctr+16) = i+ii + col(MatStencil_j,ctr+16) = j+jj + col(MatStencil_k,ctr+16) = k+kk + col(MatStencil_c,ctr+16) = 2 + enddo; enddo; enddo + row = col + ele = ele + 1 + K_ele = 0.0 + K_ele(1 :8 ,1 :8 ) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele(9 :16,9 :16) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele(17:24,17:24) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & + materialpoint_dPdF(2,2,2,2,1,ele) + & + materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele = K_ele + & + matmul(transpose(BMatFull), & + matmul(reshape(reshape(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,ele), & + shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ + call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr) + CHKERRQ(ierr) + enddo; enddo; enddo + call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! applying boundary conditions + rows = [0, 1, 2] + diag = (C_volAvg(1,1,1,1)/delta(1)**2.0_pReal + & + C_volAvg(2,2,2,2)/delta(2)**2.0_pReal + & + C_volAvg(3,3,3,3)/delta(3)**2.0_pReal)*detJ + call MatZeroRowsColumns(Jac,size(rows),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr) + CHKERRQ(ierr) + call DMGetGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) + ele = 0 + do k = zstart, zend; do j = ystart, yend; do i = xstart, xend + ele = ele + 1 + x_scal(0:2,i,j,k) = mesh_ipCoordinates(1:3,1,ele) + enddo; enddo; enddo + call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) ! initialize to undeformed coordinates (ToDo: use ip coordinates) + call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr);CHKERRQ(ierr) ! get rigid body deformation modes + call DMRestoreGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr) + call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) + call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) + call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr) + +end subroutine grid_mech_FEM_formJacobian + + +!-------------------------------------------------------------------------------------------------- +!> @brief convergence check +!-------------------------------------------------------------------------------------------------- +subroutine grid_mech_FEM_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) +use mesh +use spectral_utilities + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_stress_tolRel, & + err_stress_tolAbs + use FEsolving, only: & + terminallyIll + + implicit none + SNES :: snes_local + PetscInt :: PETScIter + PetscReal :: & + xnorm, & ! not used + snorm, & ! not used + fnorm + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal) :: & + err_div, & + divTol, & + BCTol + +err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ + divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs) + + write(6,*) BCTol,divTol + converged: if ((totalIter >= itmin -1 .and. & + all([ err_div/divTol, & + err_BC /BCTol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then converged + reason = -1 + else converged + reason = 0 + endif converged + +!-------------------------------------------------------------------------------------------------- +! report + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) + +end subroutine grid_mech_FEM_converged + +!-------------------------------------------------------------------------------------------------- +!> @brief forwarding routine +!> @details find new boundary conditions and best F estimate for end of current timestep +!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates +!-------------------------------------------------------------------------------------------------- +subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) + use math, only: & + math_mul33x33 ,& + math_rotate_backward33 + use numerics, only: & + worldrank + use homogenization, only: & + materialpoint_F0 + use mesh, only: & + grid, & + grid3 + use CPFEM2, only: & + CPFEM_age + use spectral_utilities, only: & + utilities_updateIPcoords, & + tBoundaryCondition, & + cutBack + use FEsolving, only: & + restartWrite + + implicit none + logical, intent(in) :: & + guess + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & + stress_BC, & + deformation_BC + real(pReal), dimension(3,3), intent(in) :: & + rotation_BC + PetscErrorCode :: ierr + character(len=32) :: rankStr + PetscScalar, pointer, dimension(:,:,:,:) :: & + u_current,u_lastincrement,u_rate + + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastincrement,ierr); CHKERRQ(ierr) + + if (cutBack) then + C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? + else + !-------------------------------------------------------------------------------------------------- + ! restart information for spectral solver + !if (restartWrite) then + ! write(6,'(/,a)') ' writing converged results for restart' + ! flush(6) + + ! if (worldrank == 0_pInt) then + ! call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg)) + ! write (777,rec=1) C_volAvg; close(777) + ! call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc)) + ! write (777,rec=1) C_volAvgLastInc; close(777) + ! call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot)) + ! write (777,rec=1) F_aimDot; close(777) + ! endif + + ! write(rankStr,'(a1,i0)')'_',worldrank + ! call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file + ! write (777,rec=1) F; close (777) + ! call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file + ! write (777,rec=1) F_lastInc; close (777) + ! call IO_write_jobRealFile(777,'u_current'//trim(rankStr),size(u_current)) + ! write (777,rec=1) u_current; close (777) + ! call IO_write_jobRealFile(777,'u_lastincrement'//trim(rankStr),size(u_lastincrement)) + ! write (777,rec=1) u_lastincrement; close (777) + !endif + call CPFEM_age() ! age state and kinematics + call utilities_updateIPcoords(F) + + C_volAvgLastInc = C_volAvg + + F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aim_lastInc = F_aim + + !-------------------------------------------------------------------------------------------------- + ! calculate rate for aim + if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc) + elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * deformation_BC%values + elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + endif + + + if (guess) then + call VecWAXPY(solution_rate,-1.0,solution_lastInc,solution_current,ierr) + CHKERRQ(ierr) + call VecScale(solution_rate,1.0/timeinc_old,ierr); CHKERRQ(ierr) + else + call VecSet(solution_rate,0.0,ierr); CHKERRQ(ierr) + endif + call VecCopy(solution_current,solution_lastInc,ierr); CHKERRQ(ierr) + F_lastInc = F ! winding F forward + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + endif + +!-------------------------------------------------------------------------------------------------- +! update average and local deformation gradients + F_aim = F_aim_lastInc + F_aimDot * timeinc + call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) + + call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) + CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastincrement,ierr) + CHKERRQ(ierr) + +end subroutine grid_mech_FEM_forward + +end module grid_mech_FEM