same layout for easy diff

src/grid/grid_mech_FEM.f90

@@ -34,18 +34,15 @@ module grid_mech_FEM
       itmin, &                                                                                        !< minimum number of iterations
       itmax                                                                                           !< maximum number of iterations
     real(pReal) :: &
-      err_div, &
-      divTol, &
-      BCTol, &
       eps_div_atol, &                                                                                 !< absolute tolerance for equilibrium
       eps_div_rtol, &                                                                                 !< relative tolerance for equilibrium
       eps_stress_atol, &                                                                              !< absolute tolerance for fullfillment of stress BC
       eps_stress_rtol                                                                                 !< relative tolerance for fullfillment of stress BC
   end type tNumerics
 
-  type(tNumerics) :: num
+  type(tNumerics) :: num                                                                              ! numerics parameters. Better name?
-  logical :: &
+
-    debugRotation 
+  logical :: debugRotation
 
 !--------------------------------------------------------------------------------------------------
 ! PETSc data
@@ -72,7 +69,6 @@ module grid_mech_FEM
     P_av = 0.0_pReal                                                                                !< average 1st Piola--Kirchhoff stress
 
   character(len=:), allocatable :: incInfo                                                          !< time and increment information
-
   real(pReal), dimension(3,3,3,3) :: &
     C_volAvg = 0.0_pReal, &                                                                         !< current volume average stiffness
     C_volAvgLastInc = 0.0_pReal, &                                                                  !< previous volume average stiffness
@@ -99,7 +95,6 @@ contains
 subroutine grid_mech_FEM_init
 
   real(pReal) :: HGCoeff = 0.0e-2_pReal
-  PetscInt, dimension(0:worldsize-1) :: localK
   real(pReal), dimension(3,3) :: &
     temp33_Real = 0.0_pReal
   real(pReal), dimension(4,8) :: &
@@ -111,24 +106,25 @@ subroutine grid_mech_FEM_init
                       -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])
+  real(pReal), dimension(3,3,3,3) :: devNull
   PetscErrorCode :: ierr
+  PetscScalar, pointer, dimension(:,:,:,:) :: &
+    u_current,u_lastInc
+  PetscInt, dimension(0:worldsize-1) :: localK
   integer(HID_T) :: fileHandle, groupHandle
   character(len=pStringLen) :: &
     fileName
   class(tNode), pointer :: &
     num_grid, &
     debug_grid
-  real(pReal), dimension(3,3,3,3) :: devNull
-  PetscScalar, pointer, dimension(:,:,:,:) :: &
-  u_current,u_lastInc
 
   print'(/,a)', ' <<<+-  grid_mech_FEM init  -+>>>'; flush(6)
 
-!-----------------------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------------------------
 ! debugging options
   debug_grid => config_debug%get('grid', defaultVal=emptyList)
   debugRotation = debug_grid%contains('rotation')
- 
+
 !-------------------------------------------------------------------------------------------------
 ! read numerical parameters and do sanity checks
   num_grid => config_numerics%get('grid',defaultVal=emptyDict)
@@ -242,6 +238,7 @@ subroutine grid_mech_FEM_init
     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 restartRead
+
   materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3])                           ! set starting condition for materialpoint_stressAndItsTangent
   call utilities_updateCoords(F)
   call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, &                     ! stress field, stress avg, global average of stiffness and (min+max)/2

src/grid/grid_mech_spectral_basic.f90

@@ -96,18 +96,17 @@ subroutine grid_mech_spectral_basic_init
   real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
   real(pReal), dimension(3,3) :: &
     temp33_Real = 0.0_pReal
-  class (tNode), pointer :: &
-    num_grid, &
-    debug_grid
- 
   PetscErrorCode :: ierr
   PetscScalar, pointer, dimension(:,:,:,:) :: &
     F                                                                                               ! pointer to solution data
-  PetscInt, dimension(worldsize) :: localK
+  PetscInt, dimension(0:worldsize-1) :: localK
   integer(HID_T) :: fileHandle, groupHandle
   integer        :: fileUnit
   character(len=pStringLen) :: &
     fileName
+  class (tNode), pointer :: &
+    num_grid, &
+    debug_grid
 
   print'(/,a)', ' <<<+-  grid_mech_spectral_basic init  -+>>>'; flush(6)
 
@@ -121,7 +120,7 @@ subroutine grid_mech_spectral_basic_init
 ! debugging options
   debug_grid => config_debug%get('grid', defaultVal=emptyList)
   debugRotation = debug_grid%contains('rotation')
-  
+
 !-------------------------------------------------------------------------------------------------
 ! read numerical parameters and do sanity checks
   num_grid => config_numerics%get('grid',defaultVal=emptyDict)
@@ -140,7 +139,7 @@ subroutine grid_mech_spectral_basic_init
   if (num%eps_stress_rtol < 0.0_pReal)           call IO_error(301,ext_msg='eps_stress_rtol')
   if (num%itmax <= 1)                            call IO_error(301,ext_msg='itmax')
   if (num%itmin > num%itmax .or. num%itmin < 1)  call IO_error(301,ext_msg='itmin')
-   
+
 !--------------------------------------------------------------------------------------------------
 ! set default and user defined options for PETSc
   call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr)
@@ -157,8 +156,8 @@ subroutine grid_mech_spectral_basic_init
 ! 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            = 0
-  localK(worldrank+1) = grid3
+  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
@@ -202,8 +201,8 @@ subroutine grid_mech_spectral_basic_init
   endif restartRead
 
   materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3])                           ! set starting condition for materialpoint_stressAndItsTangent
-  call Utilities_updateCoords(reshape(F,shape(F_lastInc)))
+  call utilities_updateCoords(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
+  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
   call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                ! deassociate pointer
@@ -288,7 +287,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
   type(rotation),              intent(in) :: &
     rotation_BC
   PetscErrorCode :: ierr
-  PetscScalar, dimension(:,:,:,:), pointer :: F
+  PetscScalar, pointer, dimension(:,:,:,:) :: F
 
 !--------------------------------------------------------------------------------------------------
 ! set module wide available data
@@ -334,7 +333,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
 !--------------------------------------------------------------------------------------------------
 ! 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
               rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
   call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
 
@@ -368,7 +367,7 @@ subroutine grid_mech_spectral_basic_restartWrite
 
   call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
 
-  print'(a)', ' writing solver data required for restart to file'; flush(6)
+  print*, 'writing solver data required for restart to file'; flush(6)
 
   write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
   fileHandle  = HDF5_openFile(fileName,'w')
@@ -462,6 +461,7 @@ subroutine formResidual(in, F, &
   call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
 
   if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1                                              ! new increment
+
 !--------------------------------------------------------------------------------------------------
 ! begin of new iteration
   newIteration: if (totalIter <= PETScIter) then
@@ -484,8 +484,8 @@ subroutine formResidual(in, F, &
 
 !--------------------------------------------------------------------------------------------------
 ! stress BC handling
-  deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)
+  deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)                                         ! S = 0.0 for no bc
-  F_aim = F_aim - deltaF_aim
+  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
 
 !--------------------------------------------------------------------------------------------------
@@ -493,7 +493,7 @@ subroutine formResidual(in, F, &
   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
+  err_div = utilities_divergenceRMS()                                                               ! divRMS of tensorField_fourier for later use
   call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.))       ! convolution of Gamma and tensorField_fourier
   call utilities_FFTtensorBackward                                                                  ! FFT backward of global tensorField_fourier
 

src/grid/grid_mech_spectral_polarisation.f90

@@ -15,7 +15,6 @@ module grid_mech_spectral_polarisation
   use DAMASK_interface
   use HDF5_utilities
   use math
-  use rotations
   use spectral_utilities
   use FEsolving
   use config
@@ -108,10 +107,6 @@ subroutine grid_mech_spectral_polarisation_init
   real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
   real(pReal), dimension(3,3) :: &
     temp33_Real = 0.0_pReal
-  class (tNode), pointer :: &
-    num_grid, &
-    debug_grid
-
   PetscErrorCode :: ierr
   PetscScalar, pointer, dimension(:,:,:,:) :: &
     FandF_tau, &                                                                                    ! overall pointer to solution data
@@ -122,13 +117,16 @@ subroutine grid_mech_spectral_polarisation_init
   integer        :: fileUnit
   character(len=pStringLen) :: &
     fileName
+  class (tNode), pointer :: &
+    num_grid, &
+    debug_grid
 
   print'(/,a)', ' <<<+-  grid_mech_spectral_polarisation init  -+>>>'; flush(6)
 
   print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
   print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
 
-!------------------------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------------------------
 ! debugging options
   debug_grid => config_debug%get('grid',defaultVal=emptyList)
   debugRotation = debug_grid%contains('rotation')
@@ -136,6 +134,7 @@ subroutine grid_mech_spectral_polarisation_init
 !-------------------------------------------------------------------------------------------------
 ! read numerical parameters and do sanity checks
   num_grid => config_numerics%get('grid',defaultVal=emptyDict)
+
   num%update_gamma       = num_grid%get_asBool  ('update_gamma',    defaultVal=.false.)
   num%eps_div_atol       = num_grid%get_asFloat ('eps_div_atol',    defaultVal=1.0e-4_pReal)
   num%eps_div_rtol       = num_grid%get_asFloat ('eps_div_rtol',    defaultVal=5.0e-4_pReal)
@@ -228,8 +227,8 @@ subroutine grid_mech_spectral_polarisation_init
   endif restartRead
 
   materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3])                           ! set starting condition for materialpoint_stressAndItsTangent
-  call Utilities_updateCoords(reshape(F,shape(F_lastInc)))
+  call utilities_updateCoords(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
+  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
   call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)                        ! deassociate pointer
@@ -277,7 +276,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution)
 !--------------------------------------------------------------------------------------------------
 ! update stiffness (and gamma operator)
   S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg)
-  if (num%update_gamma) then
+  if(num%update_gamma) then
     call utilities_updateGamma(C_minMaxAvg)
     C_scale = C_minMaxAvg
     S_scale = math_invSym3333(C_minMaxAvg)
@@ -320,7 +319,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
   type(rotation), intent(in) :: &
     rotation_BC
   PetscErrorCode :: ierr
-  PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
+  PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau
   integer :: i, j, k
   real(pReal), dimension(3,3) :: F_lambda33
 
@@ -588,7 +587,7 @@ subroutine formResidual(in, FandF_tau, &
 
 !--------------------------------------------------------------------------------------------------
 ! stress BC handling
-  F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC)))                                  ! S = 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 &
                                                  -params%rotation_BC%rotate(F_av)) + &
                                  params%stress_mask  * (P_av-params%stress_BC)))                    ! mask = 0.0 for no bc
@@ -596,7 +595,7 @@ subroutine formResidual(in, FandF_tau, &
   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
+  err_div = utilities_divergenceRMS()                                                               !< root mean squared error in divergence of stress
 
 !--------------------------------------------------------------------------------------------------
 ! constructing residual
@@ -615,7 +614,7 @@ subroutine formResidual(in, FandF_tau, &
   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()
+  err_curl = utilities_curlRMS()
 
 end subroutine formResidual