diff --git a/src/grid_damage_spectral.f90 b/src/grid_damage_spectral.f90 index ba9cd31ea..0663545d3 100644 --- a/src/grid_damage_spectral.f90 +++ b/src/grid_damage_spectral.f90 @@ -1,4 +1,5 @@ !-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH !> @brief Spectral solver for nonlocal damage @@ -6,43 +7,42 @@ module grid_damage_spectral #include #include - use PETScdmda - use PETScsnes - use prec, only: & - pReal - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams - - implicit none - private - + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private !-------------------------------------------------------------------------------------------------- ! derived types - type(tSolutionParams), private :: params + type(tSolutionParams), private :: params !-------------------------------------------------------------------------------------------------- ! PETSc data - SNES, private :: damage_snes - Vec, private :: solution_vec - PetscInt, private :: xstart, xend, ystart, yend, zstart, zend - real(pReal), private, dimension(:,:,:), allocatable :: & - damage_current, & !< field of current damage - damage_lastInc, & !< field of previous damage - damage_stagInc !< field of staggered damage - + SNES, private :: damage_snes + Vec, private :: solution_vec + PetscInt, private :: xstart, xend, ystart, yend, zstart, zend + real(pReal), private, dimension(:,:,:), allocatable :: & + damage_current, & !< field of current damage + damage_lastInc, & !< field of previous damage + damage_stagInc !< field of staggered damage + !-------------------------------------------------------------------------------------------------- ! reference diffusion tensor, mobility etc. - integer, private :: totalIter = 0 !< total iteration in current increment - real(pReal), dimension(3,3), private :: D_ref - real(pReal), private :: mobility_ref - - public :: & - grid_damage_spectral_init, & - grid_damage_spectral_solution, & - grid_damage_spectral_forward - private :: & - formResidual + integer, private :: totalIter = 0 !< total iteration in current increment + real(pReal), dimension(3,3), private :: D_ref + real(pReal), private :: mobility_ref + + public :: & + grid_damage_spectral_init, & + grid_damage_spectral_solution, & + grid_damage_spectral_forward + private :: & + formResidual contains @@ -51,8 +51,6 @@ contains ! ToDo: Restart not implemented !-------------------------------------------------------------------------------------------------- subroutine grid_damage_spectral_init - use IO, only: & - IO_intOut use spectral_utilities, only: & wgt use mesh, only: & @@ -94,22 +92,21 @@ subroutine grid_damage_spectral_init 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 + 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, & - 1, 0, & !< #dof (damage phase field), ghost boundary width (domain overlap) - [grid(1)],[grid(2)],localK, & !< local grid - damage_grid,ierr) !< handle, error + 1, 0, & ! #dof (damage phase field), ghost boundary width (domain overlap) + [grid(1)],[grid(2)],localK, & ! local grid + damage_grid,ierr) ! handle, error CHKERRQ(ierr) - call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) !< connect snes to da + call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) ! connect snes to da call DMsetFromOptions(damage_grid,ierr); CHKERRQ(ierr) call DMsetUp(damage_grid,ierr); CHKERRQ(ierr) - call DMCreateGlobalVector(damage_grid,solution_vec,ierr); CHKERRQ(ierr) !< global solution vector (grid x 1, i.e. every def grad tensor) - call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,& - PETSC_NULL_SNES,ierr) !< residual vector of same shape as solution vector + call DMCreateGlobalVector(damage_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) + call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) - call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) !< pull it all together with additional CLI arguments + call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr) if (trim(snes_type) == 'vinewtonrsls' .or. & trim(snes_type) == 'vinewtonssls') then @@ -117,7 +114,7 @@ subroutine grid_damage_spectral_init call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr) call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr) - call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) !< variable bounds for variational inequalities like contact mechanics, damage etc. + call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) ! variable bounds for variational inequalities like contact mechanics, damage etc. call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) endif @@ -129,10 +126,10 @@ subroutine grid_damage_spectral_init xend = xstart + xend - 1 yend = ystart + yend - 1 zend = zstart + zend - 1 - call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr) allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) + call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! damage reference diffusion update @@ -150,7 +147,8 @@ subroutine grid_damage_spectral_init call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) end subroutine grid_damage_spectral_init - + + !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral damage scheme with internal iterations !-------------------------------------------------------------------------------------------------- @@ -170,16 +168,15 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result( timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment loadCaseTime !< remaining time of current load case - integer :: i, j, k, cell + type(tSolutionState) :: solution PetscInt ::position - PetscReal :: minDamage, maxDamage, stagNorm, solnNorm + PetscReal :: minDamage, maxDamage, stagNorm, solnNorm + PetscErrorCode :: ierr - type(tSolutionState) :: & - solution SNESConvergedReason :: reason - solution%converged =.false. + solution%converged =.false. !-------------------------------------------------------------------------------------------------- ! set module wide availabe data @@ -201,7 +198,7 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result( call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) damage_stagInc = damage_current - solution%stagConverged = stagNorm < min(err_damage_tolAbs,err_damage_tolRel*solnNorm) + solution%stagConverged = stagNorm < min(err_damage_tolAbs, err_damage_tolRel*solnNorm) !-------------------------------------------------------------------------------------------------- ! updating damage state @@ -214,8 +211,8 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result( call VecMin(solution_vec,position,minDamage,ierr); CHKERRQ(ierr) call VecMax(solution_vec,position,maxDamage,ierr); CHKERRQ(ierr) if (solution%converged) & - write(6,'(/,a)') ' ... nonlocal damage converged .....................................' - write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',& + write(6,'(/,a)') ' ... nonlocal damage converged .....................................' + write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',& minDamage, maxDamage, stagNorm write(6,'(/,a)') ' ===========================================================================' flush(6) @@ -226,55 +223,55 @@ end function grid_damage_spectral_solution !-------------------------------------------------------------------------------------------------- !> @brief spectral damage forwarding routine !-------------------------------------------------------------------------------------------------- -subroutine grid_damage_spectral_forward() - use mesh, only: & - grid, & - grid3 - use spectral_utilities, only: & - cutBack, & - wgt - use damage_nonlocal, only: & - damage_nonlocal_putNonLocalDamage, & - damage_nonlocal_getDiffusion33, & - damage_nonlocal_getMobility - - implicit none - integer :: i, j, k, cell - DM :: dm_local - PetscScalar, dimension(:,:,:), pointer :: x_scal - PetscErrorCode :: ierr +subroutine grid_damage_spectral_forward + use mesh, only: & + grid, & + grid3 + use spectral_utilities, only: & + cutBack, & + wgt + use damage_nonlocal, only: & + damage_nonlocal_putNonLocalDamage, & + damage_nonlocal_getDiffusion33, & + damage_nonlocal_getMobility + + implicit none + integer :: i, j, k, cell + DM :: dm_local + PetscScalar, dimension(:,:,:), pointer :: x_scal + PetscErrorCode :: ierr - if (cutBack) then - damage_current = damage_lastInc - damage_stagInc = damage_lastInc + if (cutBack) then + damage_current = damage_lastInc + damage_stagInc = damage_lastInc !-------------------------------------------------------------------------------------------------- ! reverting damage field state - cell = 0 - call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with - x_scal(xstart:xend,ystart:yend,zstart:zend) = damage_current - call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell) - enddo; enddo; enddo - else + cell = 0 + call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with + x_scal(xstart:xend,ystart:yend,zstart:zend) = damage_current + call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell) + enddo; enddo; enddo + else !-------------------------------------------------------------------------------------------------- ! update rate and forward last inc - damage_lastInc = damage_current - cell = 0 - D_ref = 0.0_pReal - mobility_ref = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell) - mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell) - enddo; enddo; enddo - D_ref = D_ref*wgt - call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - mobility_ref = mobility_ref*wgt - call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - endif + damage_lastInc = damage_current + cell = 0 + D_ref = 0.0_pReal + mobility_ref = 0.0_pReal + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell) + mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell) + enddo; enddo; enddo + D_ref = D_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + mobility_ref = mobility_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + endif end subroutine grid_damage_spectral_forward @@ -283,84 +280,84 @@ end subroutine grid_damage_spectral_forward !> @brief forms the spectral damage residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in,x_scal,f_scal,dummy,ierr) - use numerics, only: & - residualStiffness - use mesh, only: & - grid, & - grid3 - use math, only: & - math_mul33x3 - use spectral_utilities, only: & - scalarField_real, & - vectorField_real, & - utilities_FFTvectorForward, & - utilities_FFTvectorBackward, & - utilities_FFTscalarForward, & - utilities_FFTscalarBackward, & - utilities_fourierGreenConvolution, & - utilities_fourierScalarGradient, & - utilities_fourierVectorDivergence - use damage_nonlocal, only: & - damage_nonlocal_getSourceAndItsTangent,& - damage_nonlocal_getDiffusion33, & - damage_nonlocal_getMobility - - implicit none - DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & - in - PetscScalar, dimension( & - XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & - x_scal - PetscScalar, dimension( & - X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & - f_scal - PetscObject :: dummy - PetscErrorCode :: ierr - integer :: i, j, k, cell - real(pReal) :: phiDot, dPhiDot_dPhi, mobility - - damage_current = x_scal + use numerics, only: & + residualStiffness + use mesh, only: & + grid, & + grid3 + use math, only: & + math_mul33x3 + use spectral_utilities, only: & + scalarField_real, & + vectorField_real, & + utilities_FFTvectorForward, & + utilities_FFTvectorBackward, & + utilities_FFTscalarForward, & + utilities_FFTscalarBackward, & + utilities_fourierGreenConvolution, & + utilities_fourierScalarGradient, & + utilities_fourierVectorDivergence + use damage_nonlocal, only: & + damage_nonlocal_getSourceAndItsTangent,& + damage_nonlocal_getDiffusion33, & + damage_nonlocal_getMobility + + implicit none + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & + in + PetscScalar, dimension( & + XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & + x_scal + PetscScalar, dimension( & + X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & + f_scal + PetscObject :: dummy + PetscErrorCode :: ierr + integer :: i, j, k, cell + real(pReal) :: phiDot, dPhiDot_dPhi, mobility + + damage_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field - scalarField_real = 0.0_pReal - scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current - call utilities_FFTscalarForward() - call utilities_fourierScalarGradient() !< calculate gradient of damage field - call utilities_FFTvectorBackward() - cell = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, & - vectorField_real(1:3,i,j,k)) - enddo; enddo; enddo - call utilities_FFTvectorForward() - call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field - call utilities_FFTscalarBackward() - cell = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, damage_current(i,j,k), 1, cell) - mobility = damage_nonlocal_getMobility(1,cell) - scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & - params%timeinc*phiDot + & - mobility*damage_lastInc(i,j,k) - & - mobility*damage_current(i,j,k) + & - mobility_ref*damage_current(i,j,k) - enddo; enddo; enddo + scalarField_real = 0.0_pReal + scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current + call utilities_FFTscalarForward + call utilities_fourierScalarGradient !< calculate gradient of damage field + call utilities_FFTvectorBackward + cell = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, & + vectorField_real(1:3,i,j,k)) + enddo; enddo; enddo + call utilities_FFTvectorForward + call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field + call utilities_FFTscalarBackward + cell = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, damage_current(i,j,k), 1, cell) + mobility = damage_nonlocal_getMobility(1,cell) + scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & + params%timeinc*phiDot + & + mobility*damage_lastInc(i,j,k) - & + mobility*damage_current(i,j,k) + & + mobility_ref*damage_current(i,j,k) + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! convolution of damage field with green operator - call utilities_FFTscalarForward() - call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) - call utilities_FFTscalarBackward() - where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > damage_lastInc) & - scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_lastInc - where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < residualStiffness) & - scalarField_real(1:grid(1),1:grid(2),1:grid3) = residualStiffness + call utilities_FFTscalarForward + call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) + call utilities_FFTscalarBackward + where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > damage_lastInc) & + scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_lastInc + where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < residualStiffness) & + scalarField_real(1:grid(1),1:grid(2),1:grid3) = residualStiffness !-------------------------------------------------------------------------------------------------- ! constructing residual - f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - damage_current + f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - damage_current end subroutine formResidual diff --git a/src/grid_mech_spectral_basic.f90 b/src/grid_mech_spectral_basic.f90 index 0bd04b08a..ebcc28b5e 100644 --- a/src/grid_mech_spectral_basic.f90 +++ b/src/grid_mech_spectral_basic.f90 @@ -7,65 +7,64 @@ module grid_mech_spectral_basic #include #include - use PETScdmda - use PETScsnes - use prec, only: & - pInt, & - pReal - use math, only: & - math_I3 - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams - - implicit none - private - - character (len=*), parameter, public :: & - GRID_MECH_SPECTRAL_BASIC_LABEL = 'basic' + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private + + character (len=*), parameter, public :: & + GRID_MECH_SPECTRAL_BASIC_LABEL = 'basic' !-------------------------------------------------------------------------------------------------- ! derived types - type(tSolutionParams), private :: params + type(tSolutionParams), private :: params !-------------------------------------------------------------------------------------------------- ! 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, Fdot + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot !-------------------------------------------------------------------------------------------------- ! 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 - 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) - - real(pReal), private :: & - err_BC, & !< deviation from stress BC - err_div !< RMS of div of P - - integer(pInt), private :: & - totalIter = 0_pInt !< total iteration in current increment - - public :: & - grid_mech_spectral_basic_init, & - grid_mech_spectral_basic_solution, & - grid_mech_spectral_basic_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_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 + 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) + + real(pReal), private :: & + err_BC, & !< deviation from stress BC + err_div !< RMS of div of P + + integer, private :: & + totalIter = 0 !< total iteration in current increment + + public :: & + grid_mech_spectral_basic_init, & + grid_mech_spectral_basic_solution, & + grid_mech_spectral_basic_forward contains @@ -73,145 +72,145 @@ contains !> @brief allocates all necessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine grid_mech_spectral_basic_init - use IO, only: & - IO_intOut, & - IO_error, & - IO_open_jobFile_binary - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRestart - use FEsolving, only: & - restartInc - use numerics, only: & - worldrank, & - worldsize, & - petsc_options - use homogenization, only: & - materialpoint_F0 - use DAMASK_interface, only: & - getSolverJobName - use spectral_utilities, only: & - utilities_constitutiveResponse, & - utilities_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(:,:,:,:) :: F - PetscInt, dimension(worldsize) :: localK - integer :: fileUnit - character(len=1024) :: rankStr - - write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>' - - write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' - write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012' - - 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' + use IO, only: & + IO_intOut, & + IO_error, & + IO_open_jobFile_binary + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRestart + use FEsolving, only: & + restartInc + use numerics, only: & + worldrank, & + worldsize, & + petsc_options + use homogenization, only: & + materialpoint_F0 + use DAMASK_interface, only: & + getSolverJobName + use spectral_utilities, only: & + utilities_constitutiveResponse, & + utilities_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(:,:,:,:) :: F + PetscInt, dimension(worldsize) :: localK + integer :: fileUnit + character(len=1024) :: rankStr + + write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>' + + write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' + write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012' + + 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' !-------------------------------------------------------------------------------------------------- ! 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_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) !-------------------------------------------------------------------------------------------------- ! 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, & - 9, 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 9, i.e. every def grad tensor) - call DMDASNESsetFunctionLocal(da,INSERT_VALUES,grid_mech_spectral_basic_formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector - CHKERRQ(ierr) - call SNESsetConvergenceTest(snes,grid_mech_spectral_basic_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, & + 9, 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 9, i.e. every def grad tensor) + call DMDASNESsetFunctionLocal(da,INSERT_VALUES,grid_mech_spectral_basic_formResidual,PETSC_NULL_SNES,ierr)! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESsetConvergenceTest(snes,grid_mech_spectral_basic_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,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with - - restart: if (restartInc > 0_pInt) then - if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & - 'reading values of increment ', restartInc, ' from file' - flush(6) - endif - - 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) - - F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F - call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim') - F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc - call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim_lastInc') - elseif (restartInc == 0_pInt) then restart - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = reshape(F_lastInc,[9,grid(1),grid(2),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(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,F,ierr); CHKERRQ(ierr) ! write data back to PETSc + call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with + + restart: if (restartInc > 0) then + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading values of increment ', restartInc, ' from file' + flush(6) + endif + + 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) + + F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F + call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0) call IO_error(894, ext_msg='F_aim') + F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc + call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0) call IO_error(894, ext_msg='F_aim_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]) + 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,F,ierr); CHKERRQ(ierr) ! write data back to PETSc ! QUESTION: why not writing back right after reading (l.189)? - - restartRead: if (restartInc > 0_pInt) then - if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) & - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & - 'reading more values of increment ', restartInc, ' from file' - flush(6) - 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 + + restartRead: if (restartInc > 0) then + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) & + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading more values of increment ', restartInc, ' from file' + flush(6) + 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.) @@ -222,63 +221,63 @@ end subroutine grid_mech_spectral_basic_init !> @brief solution for the Basic scheme with internal iterations !-------------------------------------------------------------------------------------------------- function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) - use numerics, only: & - update_gamma - use spectral_utilities, only: & - tBoundaryCondition, & - utilities_maskedCompliance, & - utilities_updateGamma - use FEsolving, only: & - restartWrite, & - terminallyIll - - implicit none + use numerics, only: & + update_gamma + use spectral_utilities, only: & + tBoundaryCondition, & + utilities_maskedCompliance, & + utilities_updateGamma + use FEsolving, only: & + restartWrite, & + terminallyIll + + implicit none !-------------------------------------------------------------------------------------------------- ! input data for solution - character(len=*), intent(in) :: & - incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC - type(tSolutionState) :: & - solution + 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 (update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite) + S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + if (update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite) !-------------------------------------------------------------------------------------------------- ! 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_basic_solution @@ -287,93 +286,90 @@ end function grid_mech_spectral_basic_solution !-------------------------------------------------------------------------------------------------- !> @brief forms the basic residual vector !-------------------------------------------------------------------------------------------------- -subroutine grid_mech_spectral_basic_formResidual(in, & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work) - F, & ! defgrad field on grid - residuum, & ! residuum field on grid - dummy, & - ierr) - use numerics, only: & - itmax, & - itmin - use mesh, only: & - grid, & - grid3 - use math, only: & - math_rotate_backward33, & - math_mul3333xx33 - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation - use spectral_utilities, only: & - tensorField_real, & - utilities_FFTtensorForward, & - utilities_fourierGammaConvolution, & - utilities_FFTtensorBackward, & - utilities_constitutiveResponse, & - utilities_divergenceRMS - use IO, only: & - IO_intOut - use FEsolving, only: & - terminallyIll +subroutine grid_mech_spectral_basic_formResidual(in, F, & + residuum, dummy, ierr) + use numerics, only: & + itmax, & + itmin + use mesh, only: & + grid, & + grid3 + use math, only: & + math_rotate_backward33, & + math_mul3333xx33 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation + use spectral_utilities, only: & + tensorField_real, & + utilities_FFTtensorForward, & + utilities_fourierGammaConvolution, & + utilities_FFTtensorBackward, & + utilities_constitutiveResponse, & + utilities_divergenceRMS + use IO, only: & + IO_intOut + use FEsolving, only: & + terminallyIll - implicit none - DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in - PetscScalar, & - dimension(3,3, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: F - PetscScalar, & - dimension(3,3, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: residuum - real(pReal), dimension(3,3) :: & - deltaF_aim - PetscInt :: & - PETScIter, & - nfuncs - PetscObject :: dummy - PetscErrorCode :: ierr + 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 + PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), & + intent(out) :: residuum !< residuum field + real(pReal), dimension(3,3) :: & + deltaF_aim + PetscInt :: & + PETScIter, & + nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(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, '≤', 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 !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - 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) + 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) !-------------------------------------------------------------------------------------------------- ! stress BC handling - 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 + 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 !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform - call utilities_FFTtensorForward() ! FFT forward of global "tensorField_real" - err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use - call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg - call utilities_FFTtensorBackward() ! FFT backward of global tensorField_fourier + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform + call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" + err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use + call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg + call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier !-------------------------------------------------------------------------------------------------- ! constructing residual - residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too + residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too end subroutine grid_mech_spectral_basic_formResidual @@ -382,53 +378,53 @@ end subroutine grid_mech_spectral_basic_formResidual !> @brief convergence check !-------------------------------------------------------------------------------------------------- subroutine grid_mech_spectral_basic_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) - use numerics, only: & - itmax, & - itmin, & - err_div_tolRel, & - err_div_tolAbs, & - err_stress_tolRel, & - err_stress_tolAbs - use FEsolving, only: & - terminallyIll + 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 ! not used - SNESConvergedReason :: reason - PetscObject :: dummy - PetscErrorCode :: ierr - real(pReal) :: & - divTol, & - BCTol + 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) :: & + divTol, & + BCTol - divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs) - BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_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) - converged: if ((totalIter >= itmin .and. & - all([ err_div/divTol, & - err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then - reason = 1 - elseif (totalIter >= itmax) then converged - reason = -1 - else converged - reason = 0 - endif converged + converged: if ((totalIter >= itmin .and. & + all([ err_div/divTol, & + err_BC /BCTol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then converged + reason = -1 + else converged + reason = 0 + endif converged !-------------------------------------------------------------------------------------------------- ! report - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' - write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & - err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) end subroutine grid_mech_spectral_basic_converged @@ -507,7 +503,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi write(fileUnit) F_lastInc; close (fileUnit) endif - call CPFEM_age() ! age state and kinematics + call CPFEM_age ! age state and kinematics call utilities_updateIPcoords(F) C_volAvgLastInc = C_volAvg diff --git a/src/grid_thermal_spectral.f90 b/src/grid_thermal_spectral.f90 index d1fdcdb3f..69e23c86b 100644 --- a/src/grid_thermal_spectral.f90 +++ b/src/grid_thermal_spectral.f90 @@ -1,4 +1,5 @@ !-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH !> @brief Spectral solver for thermal conduction @@ -6,43 +7,42 @@ module grid_thermal_spectral #include #include - use PETScdmda - use PETScsnes - use prec, only: & - pReal - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams - - implicit none - private - + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private !-------------------------------------------------------------------------------------------------- ! derived types - type(tSolutionParams), private :: params + type(tSolutionParams), private :: params !-------------------------------------------------------------------------------------------------- ! PETSc data - SNES, private :: thermal_snes - Vec, private :: solution_vec - PetscInt, private :: xstart, xend, ystart, yend, zstart, zend - real(pReal), private, dimension(:,:,:), allocatable :: & - temperature_current, & !< field of current temperature - temperature_lastInc, & !< field of previous temperature - temperature_stagInc !< field of staggered temperature + SNES, private :: thermal_snes + Vec, private :: solution_vec + PetscInt, private :: xstart, xend, ystart, yend, zstart, zend + real(pReal), private, dimension(:,:,:), allocatable :: & + temperature_current, & !< field of current temperature + temperature_lastInc, & !< field of previous temperature + temperature_stagInc !< field of staggered temperature !-------------------------------------------------------------------------------------------------- ! reference diffusion tensor, mobility etc. - integer, private :: totalIter = 0 !< total iteration in current increment - real(pReal), dimension(3,3), private :: D_ref - real(pReal), private :: mobility_ref - - public :: & - grid_thermal_spectral_init, & - grid_thermal_spectral_solution, & - grid_thermal_spectral_forward - private :: & - formResidual + integer, private :: totalIter = 0 !< total iteration in current increment + real(pReal), dimension(3,3), private :: D_ref + real(pReal), private :: mobility_ref + + public :: & + grid_thermal_spectral_init, & + grid_thermal_spectral_solution, & + grid_thermal_spectral_forward + private :: & + formResidual contains @@ -70,7 +70,7 @@ subroutine grid_thermal_spectral_init petsc_options implicit none - integer, dimension(worldsize) :: localK + PetscInt, dimension(worldsize) :: localK integer :: i, j, k, cell DM :: thermal_grid PetscScalar, dimension(:,:,:), pointer :: x_scal @@ -100,9 +100,9 @@ subroutine grid_thermal_spectral_init DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point grid(1),grid(2),grid(3), & ! global grid 1, 1, worldsize, & - 1, 0, & !< #dof (thermal phase field), ghost boundary width (domain overlap) - [grid(1)],[grid(2)],localK, & !< local grid - thermal_grid,ierr) !< handle, error + 1, 0, & ! #dof (thermal phase field), ghost boundary width (domain overlap) + [grid(1)],[grid(2)],localK, & ! local grid + thermal_grid,ierr) ! handle, error CHKERRQ(ierr) call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da call DMsetFromOptions(thermal_grid,ierr); CHKERRQ(ierr) @@ -123,7 +123,7 @@ subroutine grid_thermal_spectral_init allocate(temperature_lastInc(grid(1),grid(2),grid3), source=0.0_pReal) allocate(temperature_stagInc(grid(1),grid(2),grid3), source=0.0_pReal) cell = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 temperature_current(i,j,k) = temperature(material_homogenizationAt(cell))% & p(thermalMapping(material_homogenizationAt(cell))%p(1,cell)) @@ -151,6 +151,7 @@ subroutine grid_thermal_spectral_init call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) end subroutine grid_thermal_spectral_init + !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral thermal scheme with internal iterations @@ -167,7 +168,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old,loadCaseTime) result thermal_conduction_putTemperatureAndItsRate implicit none - real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment @@ -229,61 +229,61 @@ end function grid_thermal_spectral_solution !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine !-------------------------------------------------------------------------------------------------- -subroutine grid_thermal_spectral_forward() - use mesh, only: & - grid, & - grid3 - use spectral_utilities, only: & - cutBack, & - wgt - use thermal_conduction, only: & - thermal_conduction_putTemperatureAndItsRate, & - thermal_conduction_getConductivity33, & - thermal_conduction_getMassDensity, & - thermal_conduction_getSpecificHeat - - implicit none - integer :: i, j, k, cell - DM :: dm_local - PetscScalar, dimension(:,:,:), pointer :: x_scal - PetscErrorCode :: ierr - - if (cutBack) then - temperature_current = temperature_lastInc - temperature_stagInc = temperature_lastInc +subroutine grid_thermal_spectral_forward + use mesh, only: & + grid, & + grid3 + use spectral_utilities, only: & + cutBack, & + wgt + use thermal_conduction, only: & + thermal_conduction_putTemperatureAndItsRate, & + thermal_conduction_getConductivity33, & + thermal_conduction_getMassDensity, & + thermal_conduction_getSpecificHeat + + implicit none + integer :: i, j, k, cell + DM :: dm_local + PetscScalar, dimension(:,:,:), pointer :: x_scal + PetscErrorCode :: ierr + + if (cutBack) then + temperature_current = temperature_lastInc + temperature_stagInc = temperature_lastInc !-------------------------------------------------------------------------------------------------- ! reverting thermal field state - cell = 0 - call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with - x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current - call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), & - (temperature_current(i,j,k) - & - temperature_lastInc(i,j,k))/params%timeinc, & - 1,cell) - enddo; enddo; enddo - else + cell = 0 + call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with + x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current + call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), & + (temperature_current(i,j,k) - & + temperature_lastInc(i,j,k))/params%timeinc, & + 1,cell) + enddo; enddo; enddo + else !-------------------------------------------------------------------------------------------------- ! update rate and forward last inc - temperature_lastInc = temperature_current - cell = 0 - D_ref = 0.0_pReal - mobility_ref = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - D_ref = D_ref + thermal_conduction_getConductivity33(1,cell) - mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* & - thermal_conduction_getSpecificHeat(1,cell) - enddo; enddo; enddo - D_ref = D_ref*wgt - call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - mobility_ref = mobility_ref*wgt - call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - endif + temperature_lastInc = temperature_current + cell = 0 + D_ref = 0.0_pReal + mobility_ref = 0.0_pReal + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + D_ref = D_ref + thermal_conduction_getConductivity33(1,cell) + mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* & + thermal_conduction_getSpecificHeat(1,cell) + enddo; enddo; enddo + D_ref = D_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + mobility_ref = mobility_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + endif end subroutine grid_thermal_spectral_forward @@ -292,79 +292,79 @@ end subroutine grid_thermal_spectral_forward !> @brief forms the spectral thermal residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in,x_scal,f_scal,dummy,ierr) - use mesh, only: & - grid, & - grid3 - use math, only: & - math_mul33x3 - use spectral_utilities, only: & - scalarField_real, & - vectorField_real, & - utilities_FFTvectorForward, & - utilities_FFTvectorBackward, & - utilities_FFTscalarForward, & - utilities_FFTscalarBackward, & - utilities_fourierGreenConvolution, & - utilities_fourierScalarGradient, & - utilities_fourierVectorDivergence - use thermal_conduction, only: & - thermal_conduction_getSourceAndItsTangent, & - thermal_conduction_getConductivity33, & - thermal_conduction_getMassDensity, & - thermal_conduction_getSpecificHeat + use mesh, only: & + grid, & + grid3 + use math, only: & + math_mul33x3 + use spectral_utilities, only: & + scalarField_real, & + vectorField_real, & + utilities_FFTvectorForward, & + utilities_FFTvectorBackward, & + utilities_FFTscalarForward, & + utilities_FFTscalarBackward, & + utilities_fourierGreenConvolution, & + utilities_fourierScalarGradient, & + utilities_fourierVectorDivergence + use thermal_conduction, only: & + thermal_conduction_getSourceAndItsTangent, & + thermal_conduction_getConductivity33, & + thermal_conduction_getMassDensity, & + thermal_conduction_getSpecificHeat + + implicit none + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & + in + PetscScalar, dimension( & + XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & + x_scal + PetscScalar, dimension( & + X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & + f_scal + PetscObject :: dummy + PetscErrorCode :: ierr + integer :: i, j, k, cell + real(pReal) :: Tdot, dTdot_dT - implicit none - DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & - in - PetscScalar, dimension( & - XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & - x_scal - PetscScalar, dimension( & - X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & - f_scal - PetscObject :: dummy - PetscErrorCode :: ierr - integer :: i, j, k, cell - real(pReal) :: Tdot, dTdot_dT - - temperature_current = x_scal + temperature_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field - scalarField_real = 0.0_pReal - scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current - call utilities_FFTscalarForward() - call utilities_fourierScalarGradient() !< calculate gradient of damage field - call utilities_FFTvectorBackward() - cell = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, & - vectorField_real(1:3,i,j,k)) - enddo; enddo; enddo - call utilities_FFTvectorForward() - call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field - call utilities_FFTscalarBackward() - cell = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell) - scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & - params%timeinc*Tdot + & - thermal_conduction_getMassDensity (1,cell)* & - thermal_conduction_getSpecificHeat(1,cell)*(temperature_lastInc(i,j,k) - & - temperature_current(i,j,k)) + & - mobility_ref*temperature_current(i,j,k) - enddo; enddo; enddo + scalarField_real = 0.0_pReal + scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current + call utilities_FFTscalarForward + call utilities_fourierScalarGradient !< calculate gradient of damage field + call utilities_FFTvectorBackward + cell = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, & + vectorField_real(1:3,i,j,k)) + enddo; enddo; enddo + call utilities_FFTvectorForward + call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field + call utilities_FFTscalarBackward + cell = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell) + scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & + params%timeinc*Tdot + & + thermal_conduction_getMassDensity (1,cell)* & + thermal_conduction_getSpecificHeat(1,cell)*(temperature_lastInc(i,j,k) - & + temperature_current(i,j,k)) + & + mobility_ref*temperature_current(i,j,k) + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! convolution of damage field with green operator - call utilities_FFTscalarForward() - call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) - call utilities_FFTscalarBackward() + call utilities_FFTscalarForward + call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) + call utilities_FFTscalarBackward !-------------------------------------------------------------------------------------------------- ! constructing residual - f_scal = temperature_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) + f_scal = temperature_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) end subroutine formResidual