diff --git a/VERSION b/VERSION index 4d6205d89..9096149dc 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-138-g7b04b761 +v2.0.3-198-g1c762860 diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 45e423fe3..0ecea2b44 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -257,7 +257,7 @@ subroutine CPFEM_age() write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file' write(rankStr,'(a1,i0)')'_',worldrank - fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a') call HDF5_write(fileHandle,material_phase, 'recordedPhase') call HDF5_write(fileHandle,crystallite_F0, 'convergedF') diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index c5e9a254b..3ce37c5ff 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -64,7 +64,6 @@ subroutine grid_damage_spectral_init worldsize, & petsc_options - implicit none PetscInt, dimension(worldsize) :: localK integer :: i, j, k, cell DM :: damage_grid @@ -164,7 +163,6 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result( use damage_nonlocal, only: & damage_nonlocal_putNonLocalDamage - implicit none real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment @@ -236,7 +234,6 @@ subroutine grid_damage_spectral_forward damage_nonlocal_getDiffusion33, & damage_nonlocal_getMobility - implicit none integer :: i, j, k, cell DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal @@ -301,7 +298,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) damage_nonlocal_getDiffusion33, & damage_nonlocal_getMobility - implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( & diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index e31d93637..2feec76df 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -7,65 +7,66 @@ 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 + use DAMASK_interface + use HDF5_utilities + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private !-------------------------------------------------------------------------------------------------- ! derived types - type(tSolutionParams), private :: params + type(tSolutionParams), private :: params !-------------------------------------------------------------------------------------------------- ! PETSc data - DM, private :: mech_grid - SNES, private :: mech_snes - Vec, private :: solution_current, solution_lastInc, solution_rate + 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 + 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 + 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, private :: & + totalIter = 0 !< total iteration in current increment + + public :: & + grid_mech_FEM_init, & + grid_mech_FEM_solution, & + grid_mech_FEM_forward contains @@ -77,172 +78,163 @@ subroutine grid_mech_FEM_init IO_intOut, & IO_error, & IO_open_jobFile_binary - 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 - integer :: fileUnit - character(len=1024) :: rankStr - real(pReal), dimension(3,3,3,3) :: devNull - PetscScalar, pointer, dimension(:,:,:,:) :: & - u_current,u_lastInc - - write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>' + 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 + + 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 :: rank + integer(HID_T) :: fileHandle + character(len=1024) :: rankStr + real(pReal), dimension(3,3,3,3) :: devNull + PetscScalar, pointer, dimension(:,:,:,:) :: & + u_current,u_lastInc + + 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) + 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) + 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_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,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,formResidual,PETSC_NULL_SNES,ierr) - CHKERRQ(ierr) - call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr) - CHKERRQ(ierr) - call SNESSetConvergenceTest(mech_snes,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 + 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_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,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,formResidual,PETSC_NULL_SNES,ierr) + CHKERRQ(ierr) + call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr) + CHKERRQ(ierr) + call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" + CHKERRQ(ierr) + 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_pReal,ierr);CHKERRQ(ierr) - call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr) - call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,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 + call VecSet(solution_current,0.0_pReal,ierr);CHKERRQ(ierr) + call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr) + call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,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 !-------------------------------------------------------------------------------------------------- ! init fields - restart: if (restartInc > 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading values of increment ', restartInc, ' from file' - - fileUnit = IO_open_jobFile_binary('F_aim') - read(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc') - read(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot') - read(fileUnit) F_aimDot; close(fileUnit) - - write(rankStr,'(a1,i0)')'_',worldrank - - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) - read(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) - read(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u'//trim(rankStr)) - read(fileUnit) u_current; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u_lastInc'//trim(rankStr)) - read(fileUnit) u_lastInc; close (fileUnit) - - elseif (restartInc == 0) 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_lastInc,ierr) - CHKERRQ(ierr) - - restartRead: if (restartInc > 0_pInt) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('C_volAvg') - read(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - read(fileUnit) C_volAvgLastInc; close(fileUnit) - endif restartRead + restart: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading values of increment ', restartInc, ' from file' + + write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') + + call HDF5_read(fileHandle,F_aim, 'F_aim') + call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_read(fileHandle,F, 'F') + call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_read(fileHandle,u_current, 'u') + call HDF5_read(fileHandle,u_lastInc, 'u_lastInc') + + elseif (restartInc == 0) 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_lastInc,ierr) + CHKERRQ(ierr) + + restartRead: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_closeFile(fileHandle) + endif restartRead end subroutine grid_mech_FEM_init @@ -251,60 +243,57 @@ 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 + use IO, only: & + IO_error + use spectral_utilities, only: & + tBoundaryCondition, & + utilities_maskedCompliance + use FEsolving, only: & + restartWrite, & + terminallyIll !-------------------------------------------------------------------------------------------------- ! 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 - + 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 + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + 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 + 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) + 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. + call SNESGetConvergedReason(mech_snes,reason,ierr);CHKERRQ(ierr) + + solution%converged = reason > 0 + solution%iterationsNeeded = totalIter + solution%termIll = terminallyIll + terminallyIll = .false. end function grid_mech_FEM_solution @@ -335,61 +324,52 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat use FEsolving, only: & restartWrite - implicit none logical, intent(in) :: & guess real(pReal), intent(in) :: & timeinc_old, & timeinc, & - loadCaseTime !< remaining time of current load case + 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 - integer :: fileUnit + integer(HID_T) :: fileHandle character(len=32) :: rankStr PetscScalar, pointer, dimension(:,:,:,:) :: & u_current,u_lastInc - call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) if (cutBack) then C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? else !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver - if (restartWrite) then ! QUESTION: where is this logical properly set? - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - - if (worldrank == 0) then - fileUnit = IO_open_jobFile_binary('C_volAvg','w') - write(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - write(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim','w') - write(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') - write(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot','w') - write(fileUnit) F_aimDot; close(fileUnit) - endif - + if (restartWrite) then + write(6,'(/,a)') ' writing converged results for restart';flush(6) + write(rankStr,'(a1,i0)')'_',worldrank - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') - write(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') - write(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u'//trim(rankStr),'w') - write(fileUnit) u_current; close (fileUnit) - fileUnit = IO_open_jobFile_binary('u_lastInc'//trim(rankStr),'w') - write(fileUnit) u_lastInc; close (fileUnit) + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F, 'F') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_write(fileHandle,u_current, 'u') + call HDF5_write(fileHandle,u_lastInc, 'u_lastInc') + + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + + call HDF5_closeFile(fileHandle) endif - call CPFEM_age() ! age state and kinematics + call CPFEM_age ! age state and kinematics call utilities_updateIPcoords(F) C_volAvgLastInc = C_volAvg @@ -399,13 +379,13 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat !-------------------------------------------------------------------------------------------------- ! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F F_aimDot = & F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) - elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + 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 + 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 @@ -420,20 +400,20 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat 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 + 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) + 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_lastInc,ierr) - CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) + CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) + CHKERRQ(ierr) end subroutine grid_mech_FEM_forward @@ -441,59 +421,57 @@ end subroutine grid_mech_FEM_forward !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- -subroutine 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 +subroutine converged(snes_local,PETScIter,devNull1,devNull2,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 + SNES :: snes_local + PetscInt, intent(in) :: PETScIter + PetscReal, intent(in) :: & + devNull1, & + devNull2, & + 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) + 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) - - if ((totalIter >= itmin .and. & - all([ err_div/divTol, & - err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then - reason = 1 - elseif (totalIter >= itmax) then - reason = -1 - else - reason = 0 - endif + if ((totalIter >= itmin .and. & + all([ err_div/divTol, & + err_BC /BCTol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then + reason = -1 + else + reason = 0 + endif !-------------------------------------------------------------------------------------------------- ! 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) + 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 converged @@ -501,136 +479,136 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- !> @brief forms the residual vector !-------------------------------------------------------------------------------------------------- -subroutine 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 +subroutine 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 + 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) + 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 + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! 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+1, '≤', 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 + newIteration: if (totalIter <= PETScIter) then + totalIter = totalIter + 1 + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & + trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', 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) = math_rotate_backward33(F_aim,params%rotation_BC) + transpose(matmul(BMat,x_elem)) - enddo; enddo; enddo - call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) + 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) = math_rotate_backward33(F_aim,params%rotation_BC) + 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) + 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 + 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_pReal,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) + call VecSet(f_local,0.0_pReal,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) + 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 formResidual @@ -639,97 +617,96 @@ end subroutine formResidual !> @brief forms the FEM stiffness matrix !-------------------------------------------------------------------------------------------------- subroutine 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 + use mesh, only: & + mesh_ipCoordinates + use homogenization, only: & + materialpoint_dPdF - 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) + + 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) + 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 formJacobian diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 99839e50f..2daebefbd 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -7,6 +7,8 @@ module grid_mech_spectral_basic #include #include + use DAMASK_interface + use HDF5_utilities use PETScdmda use PETScsnes use prec, only: & @@ -25,8 +27,7 @@ module grid_mech_spectral_basic type(tSolutionParams), private :: params type, private :: tNumerics - logical :: & - update_gamma !< update gamma operator with current stiffness + logical :: update_gamma !< update gamma operator with current stiffness end type tNumerics type(tNumerics) :: num ! numerics parameters. Better name? @@ -40,8 +41,8 @@ module grid_mech_spectral_basic !-------------------------------------------------------------------------------------------------- ! common pointwise data real(pReal), private, dimension(:,:,:,:,:), allocatable :: & - F_lastInc, & - Fdot + F_lastInc, & !< field of previous compatible deformation gradients + Fdot !< field of assumed rate of compatible deformation gradient !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. @@ -99,23 +100,22 @@ subroutine grid_mech_spectral_basic_init use spectral_utilities, only: & utilities_constitutiveResponse, & utilities_updateGamma, & - utilities_updateIPcoords, & - wgt + utilities_updateIPcoords use mesh, only: & grid, & grid3 use math, only: & math_invSym3333 - implicit none real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & - F ! pointer to solution data + F ! pointer to solution data PetscInt, dimension(worldsize) :: localK + integer(HID_T) :: fileHandle integer :: fileUnit character(len=1024) :: rankStr @@ -163,7 +163,7 @@ subroutine grid_mech_spectral_basic_init call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) - call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "converged" + call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" CHKERRQ(ierr) call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments @@ -174,19 +174,14 @@ subroutine grid_mech_spectral_basic_init restart: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('F_aim') - read(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc') - read(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot') - read(fileUnit) F_aimDot; close(fileUnit) - write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) - read(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) - read(fileUnit) F_lastInc; close (fileUnit) + call HDF5_read(fileHandle,F_aim, 'F_aim') + call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_read(fileHandle,F, 'F') + call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') elseif (restartInc == 0) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity @@ -203,15 +198,15 @@ subroutine grid_mech_spectral_basic_init restartRead: if (restartInc > 0) then write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('C_volAvg') - read(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - read(fileUnit) C_volAvgLastInc; close(fileUnit) + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_closeFile(fileHandle) + fileUnit = IO_open_jobFile_binary('C_ref') read(fileUnit) C_minMaxAvg; close(fileUnit) endif restartRead - call Utilities_updateGamma(C_minMaxAvg,.true.) + call utilities_updateGamma(C_minMaxAvg,.true.) end subroutine grid_mech_spectral_basic_init @@ -228,8 +223,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ restartWrite, & terminallyIll - implicit none - !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & @@ -251,8 +244,8 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) - if (num%update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite) + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg,restartWrite) !-------------------------------------------------------------------------------------------------- ! set module wide available data @@ -306,7 +299,6 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi use FEsolving, only: & restartWrite - implicit none logical, intent(in) :: & guess real(pReal), intent(in) :: & @@ -321,7 +313,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: F - integer :: fileUnit + integer(HID_T) :: fileHandle character(len=32) :: rankStr call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) @@ -331,29 +323,24 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? else !-------------------------------------------------------------------------------------------------- - ! restart information for spectral solver - if (restartWrite) then ! QUESTION: where is this logical properly set? - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - - if (worldrank == 0) then - fileUnit = IO_open_jobFile_binary('C_volAvg','w') - write(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - write(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim','w') - write(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') - write(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot','w') - write(fileUnit) F_aimDot; close(fileUnit) - endif - + ! restart information for spectral solver + if (restartWrite) then + write(6,'(/,a)') ' writing converged results for restart';flush(6) + write(rankStr,'(a1,i0)')'_',worldrank - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') - write(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') - write(fileUnit) F_lastInc; close (fileUnit) + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F, 'F') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_write(fileHandle,C_minMaxAvg, 'C_minMaxAvg') + + call HDF5_closeFile(fileHandle) endif call CPFEM_age ! age state and kinematics @@ -399,7 +386,7 @@ end subroutine grid_mech_spectral_basic_forward !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- -subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) +subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr) use numerics, only: & itmax, & itmin, & @@ -410,13 +397,12 @@ subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) use FEsolving, only: & terminallyIll - implicit none SNES :: snes_local - PetscInt :: PETScIter - PetscReal :: & - xnorm, & ! not used - snorm, & ! not used - fnorm ! not used + PetscInt, intent(in) :: PETScIter + PetscReal, intent(in) :: & + devNull1, & + devNull2, & + devNull3 SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr @@ -452,7 +438,7 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- -!> @brief forms the basic residual vector +!> @brief forms the residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in, F, & residuum, dummy, ierr) @@ -481,7 +467,6 @@ subroutine formResidual(in, F, & use FEsolving, only: & terminallyIll - implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), & intent(in) :: F !< deformation gradient field @@ -515,7 +500,7 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) + call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) P_av,C_volAvg,C_minMaxAvg, & F,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index aff4913b1..3bd30a360 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -7,78 +7,79 @@ module grid_mech_spectral_polarisation #include #include - use PETScdmda - use PETScsnes - use prec, only: & - pReal - use math, only: & - math_I3 - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams - - implicit none - private + use DAMASK_interface + use HDF5_utilities + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private !-------------------------------------------------------------------------------------------------- ! derived types type(tSolutionParams), private :: params type, private :: tNumerics - logical :: & - update_gamma !< update gamma operator with current stiffness + logical :: update_gamma !< update gamma operator with current stiffness end type tNumerics type(tNumerics) :: num ! numerics parameters. Better name? !-------------------------------------------------------------------------------------------------- ! PETSc data - DM, private :: da - SNES, private :: snes - Vec, private :: solution_vec + DM, private :: da + SNES, private :: snes + Vec, private :: solution_vec !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: & - F_lastInc, & !< field of previous compatible deformation gradients - F_tau_lastInc, & !< field of previous incompatible deformation gradient - Fdot, & !< field of assumed rate of compatible deformation gradient - F_tauDot !< field of assumed rate of incopatible deformation gradient + real(pReal), private, dimension(:,:,:,:,:), allocatable :: & + F_lastInc, & !< field of previous compatible deformation gradients + F_tau_lastInc, & !< field of previous incompatible deformation gradient + Fdot, & !< field of assumed rate of compatible deformation gradient + F_tauDot !< field of assumed rate of incopatible deformation gradient !-------------------------------------------------------------------------------------------------- ! 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_lastInc = math_I3, & !< previous average deformation gradient - F_av = 0.0_pReal, & !< average incompatible def grad field - 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 - C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness - C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness - S = 0.0_pReal, & !< current compliance (filled up with zeros) - C_scale = 0.0_pReal, & - S_scale = 0.0_pReal - - real(pReal), private :: & - err_BC, & !< deviation from stress BC - err_curl, & !< RMS of curl of F - err_div !< RMS of div of P + 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_lastInc = math_I3, & !< previous average deformation gradient + F_av = 0.0_pReal, & !< average incompatible def grad field + P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - integer, private :: & - totalIter = 0 !< total iteration in current increment + 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 + C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness + C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness + S = 0.0_pReal, & !< current compliance (filled up with zeros) + C_scale = 0.0_pReal, & + S_scale = 0.0_pReal - public :: & - grid_mech_spectral_polarisation_init, & - grid_mech_spectral_polarisation_solution, & - grid_mech_spectral_polarisation_forward - private :: & - converged, & - formResidual + real(pReal), private :: & + err_BC, & !< deviation from stress BC + err_curl, & !< RMS of curl of F + err_div !< RMS of div of P + + integer, private :: & + totalIter = 0 !< total iteration in current increment + + public :: & + grid_mech_spectral_polarisation_init, & + grid_mech_spectral_polarisation_solution, & + grid_mech_spectral_polarisation_forward + private :: & + converged, & + formResidual contains @@ -86,149 +87,141 @@ contains !> @brief allocates all necessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine grid_mech_spectral_polarisation_init - use IO, only: & - IO_intOut, & - IO_error, & - IO_open_jobFile_binary - use FEsolving, only: & - restartInc - use config, only :& - config_numerics - use numerics, only: & - worldrank, & - worldsize, & - petsc_options - use homogenization, only: & - materialpoint_F0 - use DAMASK_interface, only: & - getSolverJobName - use spectral_utilities, only: & - utilities_constitutiveResponse, & - utilities_updateGamma, & - utilities_updateIPcoords, & - wgt - use mesh, only: & - grid, & - grid3 - use math, only: & - math_invSym3333 - - implicit none - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P - real(pReal), dimension(3,3) :: & - temp33_Real = 0.0_pReal - - PetscErrorCode :: ierr - PetscScalar, pointer, dimension(:,:,:,:) :: & - FandF_tau, & ! overall pointer to solution data - F, & ! specific (sub)pointer - F_tau ! specific (sub)pointer - PetscInt, dimension(worldsize) :: localK - integer :: fileUnit - character(len=1024) :: rankStr + use IO, only: & + IO_intOut, & + IO_error, & + IO_open_jobFile_binary + use FEsolving, only: & + restartInc + use config, only :& + config_numerics + use numerics, only: & + worldrank, & + worldsize, & + petsc_options + use homogenization, only: & + materialpoint_F0 + use DAMASK_interface, only: & + getSolverJobName + use spectral_utilities, only: & + utilities_constitutiveResponse, & + utilities_updateGamma, & + utilities_updateIPcoords + use mesh, only: & + grid, & + grid3 + use math, only: & + math_invSym3333 + + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P + real(pReal), dimension(3,3) :: & + temp33_Real = 0.0_pReal - write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>' - - write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' - write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' + PetscErrorCode :: ierr + PetscScalar, pointer, dimension(:,:,:,:) :: & + FandF_tau, & ! overall pointer to solution data + F, & ! specific (sub)pointer + F_tau ! specific (sub)pointer + PetscInt, dimension(worldsize) :: localK + integer(HID_T) :: fileHandle + integer :: fileUnit + character(len=1024) :: rankStr + + write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>' - num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0 + write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' + write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' + + num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0 !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) - CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) - CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) + CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (F_tauDot (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) + allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc - call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - localK = 0 - localK(worldrank+1) = grid3 - call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) - call DMDACreate3d(PETSC_COMM_WORLD, & - DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary - DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - grid(1),grid(2),grid(3), & ! global grid - 1 , 1, worldsize, & - 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) - [grid(1)],[grid(2)],localK, & ! local grid - da,ierr) ! handle, error - CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da - call DMsetFromOptions(da,ierr); CHKERRQ(ierr) - call DMsetUp(da,ierr); CHKERRQ(ierr) - call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor) - call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector - CHKERRQ(ierr) - call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" - CHKERRQ(ierr) - call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments + call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) + localK = 0 + localK(worldrank+1) = grid3 + call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary + DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point + grid(1),grid(2),grid(3), & ! global grid + 1 , 1, worldsize, & + 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) + [grid(1)],[grid(2)],localK, & ! local grid + da,ierr) ! handle, error + CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da + call DMsetFromOptions(da,ierr); CHKERRQ(ierr) + call DMsetUp(da,ierr); CHKERRQ(ierr) + call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor) + call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" + CHKERRQ(ierr) + call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments !-------------------------------------------------------------------------------------------------- ! init fields - call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data - F => FandF_tau( 0: 8,:,:,:) - F_tau => FandF_tau( 9:17,:,:,:) + call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data + F => FandF_tau( 0: 8,:,:,:) + F_tau => FandF_tau( 9:17,:,:,:) + + restart: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file' + + write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') + + call HDF5_read(fileHandle,F_aim, 'F_aim') + call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') + call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_read(fileHandle,F, 'F') + call HDF5_read(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_read(fileHandle,F_tau, 'F_tau') + call HDF5_read(fileHandle,F_tau_lastInc,'F_tau_lastInc') + + elseif (restartInc == 0) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + F_tau = 2.0_pReal*F + F_tau_lastInc = 2.0_pReal*F_lastInc + endif restart + + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) + call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + reshape(F,shape(F_lastInc)), & ! target F + 0.0_pReal, & ! time increment + math_I3) ! no rotation of boundary condition + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer + + restartRead: if (restartInc > 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file' + call HDF5_read(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + call HDF5_closeFile(fileHandle) - restart: if (restartInc > 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file' - - fileUnit = IO_open_jobFile_binary('F_aim') - read(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc') - read(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot') - read(fileUnit) F_aimDot; close(fileUnit) - - write(rankStr,'(a1,i0)')'_',worldrank - - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) - read(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) - read(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr)) - read(fileUnit) F_tau; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr)) - read(fileUnit) F_tau_lastInc; close (fileUnit) - - elseif (restartInc == 0) then restart - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) - F_tau = 2.0_pReal*F - F_tau_lastInc = 2.0_pReal*F_lastInc - endif restart - - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 - reshape(F,shape(F_lastInc)), & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition - call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer - - restartRead: if (restartInc > 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file' - fileUnit = IO_open_jobFile_binary('C_volAvg') - read(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - read(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_ref') - read(fileUnit) C_minMaxAvg; close(fileUnit) - endif restartRead - - call Utilities_updateGamma(C_minMaxAvg,.true.) - C_scale = C_minMaxAvg - S_scale = math_invSym3333(C_minMaxAvg) + fileUnit = IO_open_jobFile_binary('C_ref') + read(fileUnit) C_minMaxAvg; close(fileUnit) + endif restartRead + + call utilities_updateGamma(C_minMaxAvg,.true.) + C_scale = C_minMaxAvg + S_scale = math_invSym3333(C_minMaxAvg) end subroutine grid_mech_spectral_polarisation_init @@ -237,66 +230,64 @@ end subroutine grid_mech_spectral_polarisation_init !> @brief solution for the Polarisation scheme with internal iterations !-------------------------------------------------------------------------------------------------- function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) - use math, only: & - math_invSym3333 - use spectral_utilities, only: & - tBoundaryCondition, & - utilities_maskedCompliance, & - utilities_updateGamma - use FEsolving, only: & - restartWrite, & - terminallyIll - - implicit none + use math, only: & + math_invSym3333 + use spectral_utilities, only: & + tBoundaryCondition, & + utilities_maskedCompliance, & + utilities_updateGamma + use FEsolving, only: & + restartWrite, & + terminallyIll !-------------------------------------------------------------------------------------------------- ! 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 + 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 + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) - if (num%update_gamma) then - call utilities_updateGamma(C_minMaxAvg,restartWrite) - C_scale = C_minMaxAvg - S_scale = math_invSym3333(C_minMaxAvg) - endif + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + if (num%update_gamma) then + call utilities_updateGamma(C_minMaxAvg,restartWrite) + C_scale = C_minMaxAvg + S_scale = math_invSym3333(C_minMaxAvg) + endif !-------------------------------------------------------------------------------------------------- ! 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 + 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(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) + call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence - call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) - - solution%converged = reason > 0 - solution%iterationsNeeded = totalIter - solution%termIll = terminallyIll - terminallyIll = .false. + call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) + + solution%converged = reason > 0 + solution%iterationsNeeded = totalIter + solution%termIll = terminallyIll + terminallyIll = .false. end function grid_mech_spectral_polarisation_solution @@ -330,13 +321,12 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa use FEsolving, only: & restartWrite - implicit none logical, intent(in) :: & guess real(pReal), intent(in) :: & timeinc_old, & timeinc, & - loadCaseTime !< remaining time of current load case + loadCaseTime !< remaining time of current load case type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC @@ -347,7 +337,7 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa integer :: i, j, k real(pReal), dimension(3,3) :: F_lambda33 - integer :: fileUnit + integer(HID_T) :: fileHandle character(len=32) :: rankStr call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) @@ -355,37 +345,29 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa F_tau => FandF_tau( 9:17,:,:,:) if (cutBack) then - C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? - C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? + C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? + C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? else !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver - if (restartWrite) then ! QUESTION: where is this logical properly set? - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - - if (worldrank == 0) then - fileUnit = IO_open_jobFile_binary('C_volAvg','w') - write(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - write(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim','w') - write(fileUnit) F_aim; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w') - write(fileUnit) F_aim_lastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot','w') - write(fileUnit) F_aimDot; close(fileUnit) - endif + if (restartWrite) then + write(6,'(/,a)') ' writing converged results for restart';flush(6) write(rankStr,'(a1,i0)')'_',worldrank - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') - write(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') - write(fileUnit) F_lastInc; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr),'w') - write(fileUnit) F_tau; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr),'w') - write(fileUnit) F_tau_lastInc; close (fileUnit) + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + + call HDF5_write(fileHandle,F_aim, 'F_aim') + call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc') + call HDF5_write(fileHandle,F_aimDot, 'F_aimDot') + call HDF5_write(fileHandle,F, 'F') + call HDF5_write(fileHandle,F_lastInc, 'F_lastInc') + call HDF5_write(fileHandle,F_tau, 'F_tau') + call HDF5_write(fileHandle,F_tau_lastInc, 'F_tau_lastInc') + + call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') + call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') + + call HDF5_closeFile(fileHandle) endif call CPFEM_age ! age state and kinematics @@ -399,13 +381,13 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa !-------------------------------------------------------------------------------------------------- ! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F F_aimDot = & F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) - elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + 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 + 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 @@ -417,33 +399,33 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa F_tauDot = utilities_calculateRate(guess, & F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, & math_rotate_backward33(F_aimDot,rotation_BC)) - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward - F_tau_lastInc = reshape(F_tau, [3,3,grid(1),grid(2),grid3]) ! winding F_tau forward - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward + F_tau_lastInc = reshape(F_tau, [3,3,grid(1),grid(2),grid3]) ! winding F_tau 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 - F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average + F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average math_rotate_backward33(F_aim,rotation_BC)),& [9,grid(1),grid(2),grid3]) - if (guess) then - F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & - [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition - else - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) - F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) - F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, & - math_mul3333xx33(C_scale,& - matmul(transpose(F_lambda33),& - F_lambda33)-math_I3))*0.5_pReal)& - + math_I3 - F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) - enddo; enddo; enddo - endif - - call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) + if (guess) then + F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & + [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition + else + do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) + F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, & + math_mul3333xx33(C_scale,& + matmul(transpose(F_lambda33),& + F_lambda33)-math_I3))*0.5_pReal)& + + math_I3 + F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) + enddo; enddo; enddo + endif + + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) end subroutine grid_mech_spectral_polarisation_forward @@ -451,208 +433,207 @@ end subroutine grid_mech_spectral_polarisation_forward !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- -subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) - use numerics, only: & - itmax, & - itmin, & - err_div_tolRel, & - err_div_tolAbs, & - err_curl_tolRel, & - err_curl_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 ! not used - SNESConvergedReason :: reason - PetscObject :: dummy - PetscErrorCode :: ierr - real(pReal) :: & - curlTol, & - divTol, & - BCTol - - curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) - divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) - BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) - - if ((totalIter >= itmin .and. & - all([ err_div /divTol, & - err_curl/curlTol, & - err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then - reason = 1 - elseif (totalIter >= itmax) then - reason = -1 - else - reason = 0 - endif +subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr) + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_curl_tolRel, & + err_curl_tolAbs, & + err_stress_tolRel, & + err_stress_tolAbs + use FEsolving, only: & + terminallyIll + + SNES :: snes_local + PetscInt, intent(in) :: PETScIter + PetscReal, intent(in) :: & + devNull1, & + devNull2, & + devNull3 + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal) :: & + curlTol, & + divTol, & + BCTol + + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) + + if ((totalIter >= itmin .and. & + all([ err_div /divTol, & + err_curl/curlTol, & + err_BC /BCTol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then + reason = -1 + else + reason = 0 + endif !-------------------------------------------------------------------------------------------------- ! 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 curl = ', & - err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' - write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & - err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) + 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 curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) end subroutine converged !-------------------------------------------------------------------------------------------------- -!> @brief forms the polarisation residual vector +!> @brief forms the residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in, FandF_tau, & residuum, dummy,ierr) - use numerics, only: & - itmax, & - itmin, & - polarAlpha, & - polarBeta - use mesh, only: & - grid, & - grid3 - use math, only: & - math_rotate_forward33, & - math_rotate_backward33, & - math_mul3333xx33, & - math_invSym3333 - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation - use spectral_utilities, only: & - wgt, & - tensorField_real, & - utilities_FFTtensorForward, & - utilities_fourierGammaConvolution, & - utilities_FFTtensorBackward, & - utilities_constitutiveResponse, & - utilities_divergenceRMS, & - utilities_curlRMS - use IO, only: & - IO_intOut - use homogenization, only: & - materialpoint_dPdF - use FEsolving, only: & - terminallyIll + use numerics, only: & + itmax, & + itmin, & + polarAlpha, & + polarBeta + use mesh, only: & + grid, & + grid3 + use math, only: & + math_rotate_forward33, & + math_rotate_backward33, & + math_mul3333xx33, & + math_invSym3333 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation + use spectral_utilities, only: & + wgt, & + tensorField_real, & + utilities_FFTtensorForward, & + utilities_fourierGammaConvolution, & + utilities_FFTtensorBackward, & + utilities_constitutiveResponse, & + utilities_divergenceRMS, & + utilities_curlRMS + use IO, only: & + IO_intOut + use homogenization, only: & + materialpoint_dPdF + use FEsolving, only: & + terminallyIll - implicit none - DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) - PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), & - target, intent(in) :: FandF_tau - PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),& - target, intent(out) :: residuum !< residuum field - PetscScalar, pointer, dimension(:,:,:,:,:) :: & - F, & - F_tau, & - residual_F, & - residual_F_tau - PetscInt :: & - PETScIter, & - nfuncs - PetscObject :: dummy - PetscErrorCode :: ierr - integer :: & - i, j, k, e + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) + PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), & + target, intent(in) :: FandF_tau + PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),& + target, intent(out) :: residuum !< residuum field + PetscScalar, pointer, dimension(:,:,:,:,:) :: & + F, & + F_tau, & + residual_F, & + residual_F_tau + PetscInt :: & + PETScIter, & + nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + integer :: & + i, j, k, e - F => FandF_tau(1:3,1:3,1,& - XG_RANGE,YG_RANGE,ZG_RANGE) - F_tau => FandF_tau(1:3,1:3,2,& - XG_RANGE,YG_RANGE,ZG_RANGE) - residual_F => residuum(1:3,1:3,1,& - X_RANGE, Y_RANGE, Z_RANGE) - residual_F_tau => residuum(1:3,1:3,2,& - X_RANGE, Y_RANGE, Z_RANGE) + F => FandF_tau(1:3,1:3,1,& + XG_RANGE,YG_RANGE,ZG_RANGE) + F_tau => FandF_tau(1:3,1:3,2,& + XG_RANGE,YG_RANGE,ZG_RANGE) + residual_F => residuum(1:3,1:3,1,& + X_RANGE, Y_RANGE, Z_RANGE) + residual_F_tau => residuum(1:3,1:3,2,& + X_RANGE, Y_RANGE, Z_RANGE) - F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt - call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt + call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) - if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment !-------------------------------------------------------------------------------------------------- ! begin of new iteration - newIteration: if (totalIter <= PETScIter) then - totalIter = totalIter + 1 - 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 + newIteration: if (totalIter <= PETScIter) then + totalIter = totalIter + 1 + 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 !-------------------------------------------------------------------------------------------------- ! - tensorField_real = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) - tensorField_real(1:3,1:3,i,j,k) = & - polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& - polarAlpha*matmul(F(1:3,1:3,i,j,k), & - math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) - enddo; enddo; enddo + tensorField_real = 0.0_pReal + do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + tensorField_real(1:3,1:3,i,j,k) = & + polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& + polarAlpha*matmul(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space - call utilities_FFTtensorForward - call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) - call utilities_FFTtensorBackward + call utilities_FFTtensorForward + call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) + residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) - P_av,C_volAvg,C_minMaxAvg, & - F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) + call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) + P_av,C_volAvg,C_minMaxAvg, & + F - residual_F_tau/polarBeta,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 = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & - -math_rotate_forward33(F_av,params%rotation_BC)) + & - params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc + err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & + -math_rotate_forward33(F_av,params%rotation_BC)) + & + params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc ! calculate divergence - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise - call utilities_FFTtensorForward - err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise + call utilities_FFTtensorForward + err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress + !-------------------------------------------------------------------------------------------------- ! constructing residual - e = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) - e = e + 1 - residual_F(1:3,1:3,i,j,k) = & - math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & - residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & - math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & - + residual_F_tau(1:3,1:3,i,j,k) - enddo; enddo; enddo + e = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + e = e + 1 + residual_F(1:3,1:3,i,j,k) = & + math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & + residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & + math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & + + residual_F_tau(1:3,1:3,i,j,k) + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! calculating curl - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F - call utilities_FFTtensorForward - err_curl = Utilities_curlRMS() + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F + call utilities_FFTtensorForward + err_curl = Utilities_curlRMS() end subroutine formResidual diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 31740feb9..e899fd89a 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -69,7 +69,6 @@ subroutine grid_thermal_spectral_init worldsize, & petsc_options - implicit none PetscInt, dimension(worldsize) :: localK integer :: i, j, k, cell DM :: thermal_grid @@ -167,7 +166,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old,loadCaseTime) result use thermal_conduction, only: & thermal_conduction_putTemperatureAndItsRate - implicit none real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment @@ -242,7 +240,6 @@ subroutine grid_thermal_spectral_forward thermal_conduction_getMassDensity, & thermal_conduction_getSpecificHeat - implicit none integer :: i, j, k, cell DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal @@ -311,7 +308,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) thermal_conduction_getMassDensity, & thermal_conduction_getSpecificHeat - implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( &