grid_damage_spectral.f90

Go to the documentation of this file.

# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/grid_damage_spectral.f90"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/grid_damage_spectral.f90"
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
module grid_damage_spectral
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscsnes.h" 1
!
!  Include file for Fortran use of the SNES package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscksp.h" 1
!
!
!  Include file for Fortran use of the KSP package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscpc.h" 1
!
!
!  Include file for Fortran use of the PC (preconditioner) package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h" 1
!
!
!  Include file for Fortran use of the Mat package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscvec.h" 1
!
!
!  Include file for Fortran use of the Vec package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscao.h" 1
!
!
!  Include file for Fortran use of the AO (application ordering) package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscis.h" 1
!
!
!  Include file for Fortran use of the IS (index set) package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 1
!
!
!  Part of the base include file for Fortran use of 1.
!  Note: This file should contain only define statements and
!  not the declaration of variables.
 
! No spaces for #defines as some compilers (PGI) also adds
! those additional spaces during preprocessing - bad for fixed format
!
 
 
 
# 1 "/opt/petsc-3.10.3/Intel-18.4-IntelMPI-2018/include/petscconf.h" 1
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
# 13 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petscversion.h" 1
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
# 17 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscviewer.h" 1
!
!  Include file for Fortran use of the PetscViewer package in 1
!
 
 
 
 
 
 
 
 
 
 
 
 
 
# 31 "/opt/petsc-3.10.3/include/petsc/finclude/petscviewer.h"
 
# 18 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscerror.h" 1
 
!
!  Include file for Fortran error codes
!    These are also in include/petscerror.h
!
 
 
 
# 24 "/opt/petsc-3.10.3/include/petsc/finclude/petscerror.h"
 
# 38 "/opt/petsc-3.10.3/include/petsc/finclude/petscerror.h"
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
# 19 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petsclog.h" 1
!
!  No includes needed for logging
# 20 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscbag.h" 1
!
!
!  Include file for Fortran use of the Bag package in 1
!
 
 
 
 
 
!
!  End of Fortran include file for the IS package in 1
 
# 21 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h" 2
 
!
! The real*8,complex*16 notatiton is used so that the
! 1 double/complex variables are not affected by
! compiler options like -r4,-r8, sometimes invoked
! by the user. NAG compiler does not like integer*4,real*8
 
# 41 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
 
 
 
# 63 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
!
 
 
 
 
 
!
# 85 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
!
 
 
 
 
 
 
!
 
 
 
 
!
 
!
 
 
!
# 128 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
# 150 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
!
!     Macro for templating between real and complex
!
# 174 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
 
 
 
!
!    Allows the matrix Fortran Kernels to work with single precision
!    matrix data structures
!
 
!
!     PetscLogDouble variables are used to contain double precision numbers
!     that are not used in the numerical computations, but rather in logging,
!     timing etc.
!
 
 
!
!     Macros for error checking
!
 
 
 
 
 
 
 
 
# 215 "/opt/petsc-3.10.3/include/petsc/finclude/petscsys.h"
 
 
 
 
 
 
 
 
 
 
 
 
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscis.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscviewer.h" 1
!
!  Include file for Fortran use of the type(tPetscViewer) package in 1
!
# 10 "/opt/petsc-3.10.3/include/petsc/finclude/petscis.h" 2
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscao.h" 2
 
 
 
 
 
 
 
 
 
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscvec.h" 2
 
 
 
 
 
 
 
 
 
 
 
 
 
# 30 "/opt/petsc-3.10.3/include/petsc/finclude/petscvec.h"
 
 
 
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h" 2
 
 
 
 
 
# 49 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h"
 
 
 
 
 
 
 
 
 
 
 
 
 
# 71 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h"
!
!  Matrix types
!
# 115 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h"
 
# 140 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h"
 
!
! MatMFFDType values
!
 
 
 
!
! MatSolverTypes
!
# 168 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h"
 
!
! GPU Storage Formats for CUSPARSE
!
 
 
 
!
! sparsity reducing ordering for STRUMPACK
!
 
 
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscpc.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h" 1
 
!
!  Include file for Fortran use of the DM package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscis.h" 1
!
!
!  Include file for Fortran use of the type(tIS) (index set) package in 1
!
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscvec.h" 1
!
!
!  Include file for Fortran use of the type(tVec) package in 1
!
# 10 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h" 1
!
!
!  Include file for Fortran use of the type(tMat) package in 1
!
# 11 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h" 2
 
 
 
 
 
 
 
 
 
 
 
# 32 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h"
 
 
# 47 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h"
 
# 10 "/opt/petsc-3.10.3/include/petsc/finclude/petscpc.h" 2
 
 
 
# 26 "/opt/petsc-3.10.3/include/petsc/finclude/petscpc.h"
!
! GAMG types
!
 
 
!
!  Various preconditioners
!
# 74 "/opt/petsc-3.10.3/include/petsc/finclude/petscpc.h"
 
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscksp.h" 2
 
 
 
 
# 21 "/opt/petsc-3.10.3/include/petsc/finclude/petscksp.h"
!
!  Various Krylov subspace methods
!
# 57 "/opt/petsc-3.10.3/include/petsc/finclude/petscksp.h"
!
!  Various Initial guesses for Krylov subspace methods
!
# 8 "/opt/petsc-3.10.3/include/petsc/finclude/petscsnes.h" 2
 
 
 
 
 
# 30 "/opt/petsc-3.10.3/include/petsc/finclude/petscsnes.h"
 
!
!  character*(80)
!
# 47 "/opt/petsc-3.10.3/include/petsc/finclude/petscsnes.h"
 
!
! character*(80)
!
 
 
 
 
 
 
 
!
! integer(kind=selected_int_kind(5))
!
 
 
 
 
 
 
!
!  character*(80)
!
# 78 "/opt/petsc-3.10.3/include/petsc/finclude/petscsnes.h"
 
 
 
# 9 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/grid_damage_spectral.f90" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h" 1
 
!
!  Include file for Fortran use of the 'da' (distributed array) package in 1
!
 
 
 
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscis.h" 1
!
!
!  Include file for Fortran use of the type(tIS) (index set) package in 1
!
# 9 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscvec.h" 1
!
!
!  Include file for Fortran use of the type(tVec) package in 1
!
# 10 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscmat.h" 1
!
!
!  Include file for Fortran use of the type(tMat) package in 1
!
# 11 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h" 2
 
# 1 "/opt/petsc-3.10.3/include/petsc/finclude/petscdm.h" 1
 
!
!  Include file for Fortran use of the type(tDM) package in 1
!
# 12 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h" 2
 
 
 
 
 
!
!   DMDA_LOCAL_INFO_SIZE is one large than the size incase the DA is larger than an integer (on 64 bit systems).
!   non-int fields are not accessiable from fortran.
!
# 39 "/opt/petsc-3.10.3/include/petsc/finclude/petscdmda.h"
 
 
 
 
 
 
 
 
 
 
 
 
 
# 10 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/grid/grid_damage_spectral.f90" 2
  use petscdmda
  use petscsnes
 
  use prec
  use spectral_utilities
  use discretization_grid
  use damage_nonlocal
  use numerics
 
  implicit none
  private
 
!--------------------------------------------------------------------------------------------------
! derived types
  type(tsolutionparams), private :: params
 
!--------------------------------------------------------------------------------------------------
! 1 data
  type(tsnes), private :: damage_snes
  type(tvec),  private :: solution_vec
  integer(kind=selected_int_kind(5)), private :: xstart, xend, ystart, yend, zstart, zend
  real(preal), private, dimension(:,:,:), allocatable :: &
    phi_current, &                                                                                  !< field of current damage
    phi_lastinc, &                                                                                  !< field of previous damage
    phi_staginc
 
!--------------------------------------------------------------------------------------------------
! reference diffusion tensor, mobility etc. 
  integer,                     private :: totaliter = 0
  real(preal), dimension(3,3), private :: k_ref
  real(preal), private                 :: mu_ref
  
  public :: &
    grid_damage_spectral_init, &
    grid_damage_spectral_solution, &
    grid_damage_spectral_forward
  private :: &
    formresidual
 
contains
 
!--------------------------------------------------------------------------------------------------
! ToDo: Restart not implemented
!--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_init
    
  integer(kind=selected_int_kind(5)), dimension(0:worldsize-1) :: localk  
  type(tdm) :: damage_grid
  type(tvec) :: ubound, lbound
  integer(kind=selected_int_kind(5)) :: ierr
  character(len=pStringLen) :: snes_type
 
  write(6,'(/,a)') ' <<<+-  grid_spectral_damage init  -+>>>'
 
  write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019'
  write(6,'(a)')   ' https://doi.org/10.1007/978-981-10-6855-3_80'
 
!--------------------------------------------------------------------------------------------------
! set default and user defined options for 1
 call petscoptionsinsertstring(petsc_null_options,'-damage_snes_type newtonls -damage_snes_mf &
                               &-damage_snes_ksp_ew -damage_ksp_type fgmres',ierr)
 if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 call petscoptionsinsertstring(petsc_null_options,trim(petsc_options),ierr)
 if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of 1
  call snescreate(petsc_comm_world,damage_snes,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call snessetoptionsprefix(damage_snes,'damage_',ierr);if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif 
  localk            = 0
  localk(worldrank) = 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, &
         1, 0, &                                                                                    ! #dof (damage phase field), ghost boundary width (domain overlap)
         [grid(1)],[grid(2)],localk, &                                                              ! local grid
         damage_grid,ierr)                                                                          ! handle, error
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call snessetdm(damage_snes,damage_grid,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif                                       ! connect snes to da
  call dmsetfromoptions(damage_grid,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmsetup(damage_grid,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call dmcreateglobalvector(damage_grid,solution_vec,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif                           ! 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
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif 
  call snessetfromoptions(damage_snes,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif                                          ! pull it all together with additional CLI arguments
  call snesgettype(damage_snes,snes_type,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  if (trim(snes_type) == 'vinewtonrsls' .or. &
      trim(snes_type) == 'vinewtonssls') then
    call dmgetglobalvector(damage_grid,lbound,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
    call dmgetglobalvector(damage_grid,ubound,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
    call vecset(lbound,0.0_preal,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
    call vecset(ubound,1.0_preal,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
    call snesvisetvariablebounds(damage_snes,lbound,ubound,ierr)                                    ! variable bounds for variational inequalities like contact mechanics, damage etc.
    call dmrestoreglobalvector(damage_grid,lbound,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
    call dmrestoreglobalvector(damage_grid,ubound,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  endif
 
!--------------------------------------------------------------------------------------------------
! init fields             
  call dmdagetcorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr)
  if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  xend = xstart + xend - 1
  yend = ystart + yend - 1
  zend = zstart + zend - 1 
  allocate(phi_current(grid(1),grid(2),grid3), source=1.0_preal)
  allocate(phi_lastinc(grid(1),grid(2),grid3), source=1.0_preal)
  allocate(phi_staginc(grid(1),grid(2),grid3), source=1.0_preal)
  call vecset(solution_vec,1.0_preal,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
!--------------------------------------------------------------------------------------------------
! damage reference diffusion update
  call updatereference
 
end subroutine grid_damage_spectral_init
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution)
 
  real(preal), intent(in) :: &
    timeinc, &                                                                                      !< increment in time for current solution
    timeinc_old
  integer :: i, j, k, cell
  type(tsolutionstate) :: solution
  integer(kind=selected_int_kind(5))  :: devnull
  real(kind=selected_real_kind(10)) :: phi_min, phi_max, stagnorm, solnnorm
  
  integer(kind=selected_int_kind(5)) :: ierr
  integer(kind=selected_int_kind(5)) :: reason
 
  solution%converged =.false.
 
!--------------------------------------------------------------------------------------------------
! set module wide availabe data 
  params%timeinc = timeinc
  params%timeincOld = timeinc_old
 
  call snessolve(damage_snes,petsc_null_vec,solution_vec,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call snesgetconvergedreason(damage_snes,reason,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
 
  if (reason < 1) then
    solution%converged = .false.
    solution%iterationsNeeded = itmax
  else
    solution%converged = .true.
    solution%iterationsNeeded = totaliter
  endif
  stagnorm = maxval(abs(phi_current - phi_staginc))
  solnnorm = maxval(abs(phi_current))
  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)
  phi_staginc = phi_current
  solution%stagConverged = stagnorm < max(err_damage_tolabs, err_damage_tolrel*solnnorm)
 
!--------------------------------------------------------------------------------------------------
! updating damage state 
  cell = 0
  do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
    cell = cell + 1
    call damage_nonlocal_putnonlocaldamage(phi_current(i,j,k),1,cell)
  enddo; enddo; enddo
 
  call vecmin(solution_vec,devnull,phi_min,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  call vecmax(solution_vec,devnull,phi_max,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
  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      = ',&
                                                        phi_min, phi_max, stagnorm
  write(6,'(/,a)') ' ==========================================================================='
  flush(6) 
 
end function grid_damage_spectral_solution
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_forward(cutBack)
   
  logical, intent(in) :: cutback
  integer                               :: i, j, k, cell
  type(tdm) :: dm_local
  real(kind=selected_real_kind(10)),  dimension(:,:,:), pointer     :: x_scal
  integer(kind=selected_int_kind(5))                              :: ierr
 
  if (cutback) then 
    phi_current = phi_lastinc
    phi_staginc = phi_lastinc
!--------------------------------------------------------------------------------------------------
! reverting damage field state 
    cell = 0
    call snesgetdm(damage_snes,dm_local,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
    call dmdavecgetarrayf90(dm_local,solution_vec,x_scal,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif                       
    x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current
    call dmdavecrestorearrayf90(dm_local,solution_vec,x_scal,ierr); if (ierr .ne. 0) then;call petscerrorf(ierr);return;endif
    do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
      cell = cell + 1                                                                           
      call damage_nonlocal_putnonlocaldamage(phi_current(i,j,k),1,cell)
    enddo; enddo; enddo
  else
    phi_lastinc = phi_current
    call updatereference
  endif  
 
end subroutine grid_damage_spectral_forward
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine formresidual(in,x_scal,f_scal,dummy,ierr)
 
  integer(kind=selected_int_kind(5)), dimension(25) :: &
    in
  real(kind=selected_real_kind(10)), dimension( &
    in(13)+1:in(13)+in(16),in(14)+1:in(14)+in(17),in(15)+1:in(15)+in(18)), intent(in) :: &
    x_scal
  real(kind=selected_real_kind(10)), dimension( &
    in(7)+1:in(7)+in(10),in(8)+1:in(8)+in(11),in(9)+1:in(9)+in(12)), intent(out) :: &
    f_scal
  integer(kind=selected_int_kind(10)) :: dummy
  integer(kind=selected_int_kind(5)) :: ierr
  integer :: i, j, k, cell
  real(preal)   :: phidot, dphidot_dphi, mobility
 
  phi_current = x_scal 
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
  scalarfield_real = 0.0_preal
  scalarfield_real(1:grid(1),1:grid(2),1:grid3) = phi_current 
  call utilities_fftscalarforward
  call utilities_fourierscalargradient                                                              
  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) = matmul(damage_nonlocal_getdiffusion(1,cell) - k_ref, &
                                         vectorfield_real(1:3,i,j,k))
  enddo; enddo; enddo
  call utilities_fftvectorforward
  call utilities_fouriervectordivergence                                                            
  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, phi_current(i,j,k), 1, cell)
    mobility = damage_nonlocal_getmobility(1,cell)
    scalarfield_real(i,j,k) = params%timeinc*(scalarfield_real(i,j,k) + phidot) &
                            + mobility*(phi_lastinc(i,j,k) - phi_current(i,j,k)) &
                            + mu_ref*phi_current(i,j,k)
  enddo; enddo; enddo
 
!--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator
  call utilities_fftscalarforward
  call utilities_fouriergreenconvolution(k_ref, mu_ref, params%timeinc)
  call utilities_fftscalarbackward
  where(scalarfield_real(1:grid(1),1:grid(2),1:grid3) > phi_lastinc) &
        scalarfield_real(1:grid(1),1:grid(2),1:grid3) = phi_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) - phi_current
 
end subroutine formresidual
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine updatereference
 
  integer :: i,j,k,cell,ierr
  
  cell = 0
  k_ref = 0.0_preal
  mu_ref = 0.0_preal
  do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
    cell = cell + 1
    k_ref  = k_ref  + damage_nonlocal_getdiffusion(1,cell)
    mu_ref = mu_ref + damage_nonlocal_getmobility(1,cell)
  enddo; enddo; enddo
  k_ref = k_ref*wgt
  call mpi_allreduce(mpi_in_place,k_ref,9,mpi_double,mpi_sum,petsc_comm_world,ierr)
  mu_ref = mu_ref*wgt
  call mpi_allreduce(mpi_in_place,mu_ref,1,mpi_double,mpi_sum,petsc_comm_world,ierr)
 
end subroutine updatereference
 
 
end module grid_damage_spectral